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SUMMARY 


This  report  describes  progress  on  a  model  that  predicts  how  Long 
wave  atmospheric  noise  is  affected  by  distortions  in  the'  earth- 
ionosphere  waveguide  caused  by  high-altitude  nuclear  bursts.  The 
model  will  also  represent  worldwide  noise  under  ambient  conditions. 

It  comprises  two  submodels:  the  source  submodel,  which  describes  the 
occurrence  rate  of  lightning  flashes  throughout  the  world;  and  the 
propagation  submodel,  which  describes  propagation  of  energy  radiated 
from  centers  of  lightning  act:.vity.  The  new  model  is  an  improvement 
over  existing  ones  in  several  respects.  First,  it  uses  recent  data  on 
the  actual  occurrence  rate  of  lightning  flashes  as  a  function  of  time 
of  day,  latitude,  and  season.  Such  data  circumvent  the  need  to  infer 
those  rates  from  records  of  thunderstorm  days.  Second,  it  includes 
extremely  low  frequency  and  low  frequency  noise.  Third,  it  uses  new 
data  on  the  altitude  and  orientation  of  individual  strokes  to  model 
transverse-electric  noise.  Finally,  the  model  incorporates  new 
propagation  algorithms  that  are  nearly  as  accurate  as  full-wave  algo¬ 
rithms,  but  they  reduce  computer  running  time  by  as  much  as  fivefold. 
It  would  be  costly  to  recalculate  mode  parameters  for  the  earth- 
ionosphere  waveguide  each  time  the  model  is  exercised.  Instead, 
attenuation  rate,  excitation  factor,  phase  velocity,  and  eigencosine 
are  precalculated  and  stored  for  a  wide  range  of  modes,  frequency, 
ground  conductivity,  and  spread-debris  nuclear  environments.  Those 
data  are  of  general  interest  and  are  presented  graphically  in  handbook 
format  in  Vol.  II. 
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PREFACE 


This  two-volume  report  describes  progress  made  during  the  first 
year  of  a  three-year  Defense  Nuclear  Agency  program  to  model  long  wave 
atmospheric  noise  in  ambient  and  nuclear  environments.  Most  of  the 
year's  efforts  have  been  spent  developing  fast-running  approxiiaate 
algorithms  to  compute  propagation  of  energy  radiated  by  lightning 
flashes.  That  portion  of  the  program  is  virtually  complete.  In 
addition,  a  literature  search  was  conducted  to  locate  recent  data  on 
the  occurrence  rate  of  worldwide  lightning  and  on  the  altitude  and 
orientation  of  typical  discharges.  Future  work  will  include  incor¬ 
porating  model  current  moments  for  several  flash  types,  combining  the 
propagation  algorithms  and  source  data  into  a  cohesive  noise  model, 
and  testing  that  model  against  available  noise  data. 
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SECTION  1 
INTRODUCTION 


Before  predictions  can  be  made  about  how  well  long-wave  communica¬ 
tion  links  would  perform  in  nuclear  environments,  it  is  necessary  to 
know  how  those  environments  affect  atmospheric  noise  as  well  as  com¬ 
munication  signals.  Because  no  data  exist  for  atmospheric  noise  under 
nuclear-disturbed  conditions,  a  predictive  noise  model  must  be 
developed  that  accounts  for  distortions  in  the  earth-ionosphere 
waveguide  caused  by  high-altitude  nuclear  bursts.  That  model  must,  of 
course,  also  accurately  represent  worldwide  noise  under  ambient 
conditions.  This  report  describes  Pacific-Sierra  Research 
Corporation's  (PSR)  progress  on  such  a  model. 

Although  man-made  or  extraterrestrial  sources  might  contribute  to 
long-wave  noise  at  certain  times  or  locations,  the  model  is  based  on 
the  assumption  that  long-wave  noise  is  caused  solely  by  lightning, 
which  radiates  strongly  in  the  extremely  low  frequency  (ELF) ,  very  low 
frequency  (VLF) ,  and  low  frequency  (LF)  bands.  The  model  therefore 
consists  of  two  submodels.  The  first,  the  source  submodel,  describes 
the  occurrence  rate  of  lightning  flashes  throughout  the  world  and  the 
electromagnetic  radiation  from  those  flashes.  The  second,  the 
propagation  submodel,  describes  how  the  energy  radiated  from  numerous 
worldwide  centers  of  lightning  activity  propagates  in  the  earth- 
ionosphere  waveguide  to  the  location  of  a  long-wave  receiver. 

Two  decades  ago,  the  Westinghouse  Georesearch  Laboratory  (WGL) 
[1972]  used  the  structure  just  described  to  develop  a  model  to  predict 
ambient  VLF  noise  throughout  the  world.  The  WGL  model  was  state  of 
the  art  at  the  time  of  its  development,  but  is  now  obsolete  in  several 
respects.  First,  data  on  actual  lightning  flash  occurrence  rates  were 
scarce,  so  WGL  inferred  the  rates  from  records  of  thunderstorm  days 
(TD) ,  thus  introducing  an  inaccuracy  of  unknown  magnitude.  Second, 

WGL  modeled  only  radiation  and  propagation  at  frequencies  between  10 
and  30  kHz,  thus  omitting  ELF  and  LF  noise  of  importance  to  systems 
deployed  during  the  past  two  decades.  Third,  the  understanding  of 
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lightning  flash  structure  has  advanced  over  the  past  20  years. 

Fourth,  the  WGL  model  neglected  horizontal  lightning  and  hence  omitted 
transverse  electric  (TE)  noise,  which  is  important  for  airborne 
receivers.  Finally,  the  propagation  algorithms  used  by  WGL  are  crude 
by  present-day  standards,  particularly  when  the  ground  or  ionosphere 
changes  along  the  direction  of  propagation;  and  they  cannot  account 
for  ionospheric  disturbances. 

Section  2  of  the  present  r.port  describes  the  results  of  a  litera¬ 
ture  search  to  locate  data  on  the  actual  occurrence  rate  of  lightning 
strokes  as  a  function  of  latitude,  time  of  day,  and  season.  Such  data 
circumvent  the  need  to  infer  lightning  occurrence  statistics  from  TDs , 
and  therefore  remove  a  major  uncertainty  from  the  modeling  process. 
Section  2  also  describes  recent  data  on  the  structure,  elevation,  and 
orientation  of  individual  strokes.  Ongoing  lightning  studies  vjill 
provide  more  data  in  the  future.  The  PSR  model  will  be  constructed  so 
such  improved  data  can  be  incorporated  easily. 

Section  3  and  Appendix  A  describe  the  propagation  algorithms  that 
will  be  used  in  the  model.  It  extends  the  work  of  Warber  and  Field 
[1987],  who  developed  approximate  methods  to  calculate  coupling  be¬ 
tween  waveguide  modes  over  geologically  complex  regions  like  Canada 
and  Greenland.  They  took  advantage  of  the  relatively  slow  lateral 
variations  of  conductivity,  and  applied  a  high-order  Wentzel-Kramers- 
Brillouin  (WKB)  approximation.  Such  approximate  methods  are  necessary 
because  full-fledged  propagation  codes  such  as  WAVEGUID  (developed  by 
the  Naval  Ocean  Systems  Center,  San  Diego,  California)  and  WAVEPROP 
(developed  by  PSR)  would  require  excessive  computer  time  if  used  for  a 
worldwide  noise  model  treating  propagation  between  hundreds  of  sources 
and  receiver  locations.  Whereas  Warber  and  Field  addressed  lateral 
transitions  in  ground  conductivity,  the  present  report  treats  changes 
in  the  waveguide  height,  such  as  between  an  undisturbed  and  nuclear- 
disturbed  region  of  the  ionosphere.  It  also  develops  a  simple  method 
to  account  for  the  anisotropy  in  propagation  caused  by  the  geomagnetic 
field,  thereby  alleviating  the  need  to  perform  separate  and  detailed 
calculations  for  each  propagation  direction  under  undisturbed 
nighttime  propagation. 
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As  an  aid,  we  have  collected  the  acronyms  used  in  the  report, 
their  definitions,  and  the  page  on  which  they  are  first  used  in  Appen¬ 
dix  B.  Appendix  C  is  a  similar  list  of  all  the  symbols  used  in  the 
mathematical  derivation  in  the  report. 

Rather  than  perform  full  calculations  each  time  the  model  is 
exercised,  it  is  preferable  to  precalculate  and  store  waveguide  mode 
parameters  for  a  wide  range  of  frequency,  ground  conductivity,  and 
spread-debris  nuclear  environments.  Those  parameters  can  be  retrieved 
rapidly,  when  needed.  The  precalculation  of  attenuation  rate,  excita¬ 
tion  factor,  phase  velocity,  and  eigen  cosine  for  all  important  modes 
has  been  performed.  Because  such  data  are  useful  in  their  own  right 
to  analysts,  they  are  presented  graphically  in  handbook  format  in 
Vol.  II  of  this  report. 
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SECTION  2 

LIGHTNING  OCCURRENCE  DATA 


Lightning  discharges  are  the  major  source  of  global  radio¬ 
frequency  (RF)  noise.  Different  types  of  lightning,  as  well  as  dif¬ 
ferent  stages  within  a  lightning  discharge,  radiate  significantly 
different  spectra.  The  noise  source  submodel  therefore  comprises 
three  major  elements:  a  representation  of  the  number  of  lightning 
discharges  for  a  given  area  during  a  specified  interval;  "modifiers," 
which  are  used  to  discriminate  between  the  various  types  of  lightning 
discharges;  and  "model  flashes,"  which  describe  the  typical  power 
spectrum  emitted  by  each  lightning  type. 

Because  of  the  large  uncertainties  inherent  in  the  WGL  data  base, 
that  model  required  a  complicated  mix  of  algorithms,  modifiers,  and  ad 
hoc  adjustments  to  produce  agreement  with  experimental  results. 
Satellite-borne  instruments  have  recently  been  used  to  measure  light¬ 
ning  occurrence  rate  and  flash  density,  which  provide  the  basis  for  a 
substantially  improved  source  submodel.  Besides  replacing  the  ob¬ 
solete  data  base,  recent  data  will  permit  a  modular  set  of  data-base 
modifiers  to  be  incorporated  in  the  new  model.  Each  modifier  will  be 
sequentially  applied  to  the  source  to  account  for  the  influence  of  an 
individual  parameter  (for  example,  storm  severity).  Although  data  to 
quantify  the  effects  of  certain  parameters — such  as  ocean  versus  land 
mass — are  still  scarce,  we  will  use  separate  modular  algorithms  that 
account  for  each  parameter,  thereby  allowing  the  modifiers  to  be 
easily  updated.  New  model  flashes  will  be  based  on  up-to-date 
theories  that  account  for  the  nonvertical ity  of  lightning  channels, 
which  is  important  since  the  new  model  must  represent  sources  of  both 
TE  and  transverse  magnetic  (TM)  noise. 

In  order  to  define  the’ terminology ,  we  first  review  lightning 
processes.  Much  of  this  review  is  taken  from  current  literature,  such 
as  Uman  and  Krider  [1982]  and  Fitzgerald  [1985].  We  then  describe  the 
old  WGL  source  maps  and  present  the  data  on  which  the  new  source  maps 
will  be  based.  Finally,  we  describe  both  the  old  and  new  source 
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modifiers.  Development  of  the  emission  spectra  for  the  various 
lightning  discharge  components  has  not  begun. 

LIGHTNING  DISCHARGE  PROCESSES. 

Lightning  is  a  complex,  propagating,  gas  breakdown  process  which 
results  in  a  transient,  high-current  electrical  discharge.  Small 
volumes  of  intense  space  charge  in  the  presence  of  a  strong  electric 
field  appear  to  initiate  the  breakdown.  A  discharge  typically 
neutralizes  some  tens  of  coulombs  of  cloud  charge  over  a  path  length 
measured  in  kilometers.  The  total  discharge,  also  called  a  "flash," 
is  made  up  of  various  components,  among  which  are  several  high-current 
pulses  called  "strokes."  A  stroke  is  defined  as  the  movement  of 
charge  in  a  given  direction.  Each  stroke  acts  as  an  individual  RF 
source,  so  basic  stroke  types  must  be  identified  and  modeled.  The 
model  flashes  are  made  up  of  combinations  of  model  strokes. 

Lightning  is  often  classified  into  one  of  two  categories.  Cloud- 
to-ground  (C-G)  discharges  (also  called  earth  flashes)  have  discharge 
channels  that  contact  the  earth.  In-cloud  (I-C)  lightning  channels  do 
not  reach  the  earth  and  are  composed  of  several  subtypes;  intercloud 
(between  two  clouds) ;  cloud-to-air  (between  a  cloud  and  clear  air) ; 
and  intracloud  (within  a  cloud) ,  which  is  by  far  the  most  common. 

Past  models  have  addressed  only  this  last  category  of  I-C  lightning, 
because  it  is  by  far  the  most  prevalent.  We  will  reevaluate  the 
assumption  that  cloud-to-air  and  intercloud  lightning  can  be  ignored. 

The  likelihood  that  a  discharge  will  produce  a  channel  to  ground 
depends  on:  the  intensity  of  the  charged  regions;  the  height  of  the 
charge  separation  aboveground;  and  the  charge  sign  which  initiates  the 
breakdown.  Other  factors,  such  as  updraft  velocity,  may  also  play  a 
role.  We  discuss  the  development  and  phases  of  I-C  and  C-G  lightning 
separately . 

In-Cloud  Lightning. 

In-cloud  lightning  flashes  occur  between  positive  and  negative 
pockets  of  space  charge.  Although  the  discharge  can  be  initiated  by 
either  charge  sign,  positive  regions  initiate  some  75  percent  of  all 
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I-C  breakdowns.  The  discharge  is  thought  to  consist  of  one  or  more 
branched,  continuously  propagating  "leaders."  The  leader  breakdowns 
generate  five  or  six  recoil  streamers  (called  "K-changes")  when  they 
contact  regions  of  space  charge  opposite  their  own.  The  development 
of  a  typical  I-C  flash  is  depicted  in  Fig.  1. 

A  typical  cloud  discharge  is  commonly  believed  to  neutralize  10 
to  30  C  of  charge  over  a  total  path  length  of  5  to  10  km  and  last 
approximately  0.5  s.  Data  from  recent  studies  of  lightning  channels 
performed  by  Richard  et  al .  [1987]  suggest  that  I-C  flashes  consist  of 
many  propagating  leaders  and  may  have  total  path  lengths  longer  than 
the  widely  accepted  5  to  10  iu...  Cloud  discharges  have  not  been 
studied  as  extensively  as  discharges  to  ground,  because  they  are  not 
as  hazardous  to  personnel  and  equipment.  Furthermore,  since  the 
charged  channels  are  at  altitude,  they  are  difficult  to  study.  As  a 
result,  less  is  known  about  them  than  C-G  discharges. 

Cloud- to-Ground  Lightning. 

The  "preliminary  breakdown"  phase  of  a  C-G  flash  initiates  in  the 
cloud  much  the  same  way  as  an  I-C  discharge,  except  that  most  C-G 
flashes  originate  from  regions  of  negative  space  charge.  The  channel 
then  extends  below  cloud  level  and  toward  ground  in  a  series  of  short 
luminous  steps  called  the  "stepped  leader."  As  the  leader  nears  the 
ground,  the  large  potential  at  the  leader  tip  induces  a  discharge 
channel  originating  from  the  ground.  The  two  channels  connect  and  the 
"main  return  stroke"  discharges  the  entire  ionized  channel.  After  the 
return-stroke  current  has  ceased  to  flow,  the  flash  may  end.  If, 
however,  "K  processes"  make  additional  charge  available  to  the  chan¬ 
nel,  a  continuous  "dart  leader"  may  propagate  down  the  primary  portion 
of  the  residual  first-stroke  channel.  The  dart  leader  differs  from 
the  stepped  leader  in  that  it  is  usually  not  branched.  The  dart 
leader  again  initiates  contact  with  the  earth,  resulting  in  another 
return  stroke.  Often  several  dart  leaders  and  subsequent  return 
strokes  occur  before  the  flash  ends.  A  slow  electric  field  change 
(called  the  "final  phase")  then  restores  the  field  in  the  vicinity  of 
the  discharge  to  ambient.  Some  leaders  begin  as  dart  leaders  but 
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a.  Step  1  — "Pre]  iminary  breakdown  leader.' 
When  leader  contacts  pocket  of  space 
charge  of  opposite  sign, 


c.  Step  3— Charge  at  source  of 
channel  causes  leader  to  propa¬ 
gate  further, 


b.  Step  2--a  "recoil  streamer"  or 
"K-change"  discharges  the  leader 
channel . 


d.  Step  4 — resulting  in  second 

K-change  which  again  discharges 
entire  leader  channel. 


Note;  Fine  lines  depict  propagating  leaders,  heavy  lines  show  discharges 

due  to  K-changss.  Notice  horizontal  structure  apparent  in  discharge 
path. 


Figure  1.  Typical  I-C  flash. 
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become  stepped  leaders  as  they  near  the  ground.  These  processes  are 
knovm  as  "dart-stepped  leaders,"  which  often  precede  the  most  intense 
return  stroke  in  a  flash.  The  phases  of  a  C-G  flash  are  shown  in 
Fig.  2. 

A  typical  C-G  discharge  lasts  between  0.1  and  1.0  s,  with  the 
majority  lasting  approximately  0.5  s.  Each  stroke  of  the  flash  lasts 
about  1  ms,  and  the  separation  time  between  strokes  is  typically  40  to 
80  ms.  The  return  stroke  generates  the  most  energetic,  natural  radio 
signals  found  on  earth,  with  a  spectral  peak  in  the  4  to  8  kHz  portion 
of  the  VLF  band.  Although  most  flashes  neutralize  negative  cloud 
charge,  positive  flashes  also  occur.  The  duration  and  intensity  of 
positive  flashes  are  usually  much  greater  than  for  negative 
discharges.  Two  model  flashes  therefore  may  be  required  to  describe 
C-G  discharges. 

NOISE  SOURCE  MAPS. 

The  noise  source  map  specifies  the  number  of  lightning  discharges 
per  unit  area  over  the  surface  of  the  earth.  Maps  are  prepared  for 
monthly  or  seasonal  periods,  and  for  specific  time  periods  over  the 
course  of  a  day.  The  well-known  CCIR  322  [1963]  set  the  standard  as 
six  4-h  daily  time  blocks  and  four  seasons.  We  next  review  the  WGL 
source  maps ,  then  present  the  data  and  format  for  the  new  global 
lightning  occurrence  (GLO)  maps. 

WGL  Thunderstorm  Day  Maps  and  Conversion  Algorithms. 

The  WGL  source  submodel  is  based  on  TD  maps  that  depict  the 
number  of  days  in  which  thunder  is  heard  in  a  particular  region.  The 
TD  data  were  provided  by  the  World  Meteorological  Organization  (WMO) . 
The  WGL  converted  other  types  of  data  into  TDs  to  improve  the  accuracy 
of  the  WMO  maps,  particularly  in  regions  where  data  were  scarce.  The 
model  also  incorporated  an  algorithm  which  attempted  to  relate  the 
number  of  TDs  in  a  month  to  a  total  number  of  lightning  discharges. 
Since  the  distribution  of  thunderstorms  over  the  course  of  a  day  is 
not  uniform,  a  table  of  hourly  diurnal  modifiers  was  then  applied  to 
allocate  the  percentages  of  a  day's  discharges  to  each  hourly  time 
block. 
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Note:  Arrows  indicate  direction  of  stroke  progression. 


The  WGL  maps  consisted  of  global  contours  of  the  number  of  TDs 
occurring  each  month.  An  empirical  relationship  (admitted  by  WGL  to 
be  imprecise)  was  proposed  to  convert  TDs  to  the  number  of  lightning 
flashes  per  square  kilometer  per  month.  That  relationship  was 


0.06  N 


1.5 

TD 


(1) 


where  =  number  of  lightning  discharges/square  kilometer/month, 

*  number  of  thunderstorm  days/month. 

The  maps  were  divided  into  sectors  of  5  deg  latitude  by  5  deg 
longitude.  The  number  of  lightning  discharges  for  each  sector  (on 
each  monthly  map)  was  calculated,  as  was  the  centroidal  location  of 
the  lightning.  Each  sector  was  then  modeled  by  a  single  noise  trans¬ 
mitter  located  at  the  lightning  centroid. 

Although  once  the  only  format  for  worldwide  data  available,  the 
TD  is  not  a  good  measure  of  lightning  activity.  Empirical  relations 
exist  for  the  average  number  of  lightning  flashes  that  occur  in  a 
typical  storm,  but  Chose  relationships  could  not  be  used  because  the 
TD  does  not  specify  Che  number  of  storms  per  day.  The  WGL  attempted 
to  correlate  the  number  of  thunderstorms  per  day  with  the  density  of 
TDs.  but  such  a  correlation  is,  at  best,  a  crude  approximation. 

New  Global  Lightning  Occurrence  Maps. 

The  advent  of  satellites  to  measure  lightning  occurrence  rates 
and  flash  density  eliminates  the  need  to  use  an  algoritbnl  to  convert 
TDs  to  lightning  occurrence.  Since  recent  satellite  studies  have 
measured  global  lightning  occurrence  at  ail  times  of  day,  diurnal 
variations  are  incorporated  into  the  raw  data,  and  no  diurnal  model  is 
required.  We  present  samples  of  such  satellite  data  and  discuss  their 
inclusion  in  the  new  source  submodel. 

OSO-2  and  OSO-5  Satellites.  The  first  satellite  studies  of 
lightning  used  optical  intensity  detectors,  which  limited  observation 
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Co  periods  of  low  light  level.  Inicial  studies,  conducted  with  Che 
OSO-2  and  OSO-5  satellites,  were  reported  by  Vorpahl,  Sparrow,  and  Ney 
[1970].  Obseir^'ations  were  restricted  to  the  2000  to  0400  time  block 
(around  local  midnight) .  and  further  limited  to  periods  of  new  moon. 
Those  studies  recorded  the  positions  of  nighttime  storms  from  February 
to  September  1969  and  January  to  July  1970.  Figure  3  shows  data  from 
OSO-5.  The  points  represent  thunderstorms,  which  must  be  converted  to 
lightning  occurrence.  This  method  is  an  improvement  over  the  TD 
because  it  provides  an  actual  storm  count.  Although  not  complete 
enough  to  fully  replace  the  TD  maps,  these  data  will  be  used  in  con¬ 
junction  with  ocher  recent  measurements,  described  below. 

Defense  Meteorological  Satellite  Program.  Subsequent  measure¬ 
ments  were  made  with  more  sophisticated  equipment,  carried  aboard  the 
Defense  Meteorological  Satellite  Program  (DMSP)  flights  2  and  3, 
Although  that  study  still  used  an  optical  intensity  detector,  it 
recorded  actual  lightning  flash  counts  at  local  dawn  and  dusk.  It 
thereby  provided  data  for  two  time  periods  and  eliminated  the  need  for 
a  conversion  between  thunderstorms  and  number  of  lightning  flashes. 
Turman  and  Edgar  [1982]  accumulated  the  data  into  10  deg  latitude  by 
10  deg  longitude  sectors  and  presented  it  as  sector  flash  rate  maps 
for  six  periods;  August,  Septeraber-October ,  November-December , 

J anuary-March ,  April,  and  Hay-June.  Figure  ka^  shows  a  sample  of  the 
flash  rate  maps  at  dawn  for  the  month  of  April. 

The  sensors  were  operated  in  a  triggered  mode  where  two  criteria 
had  to  be  met.  The  irradiance  had  to  (1)  exceed  5  x  10"^  W/cra^  (an 
unattenuated  source  power  of  4  x  10^  W)  and  (2)  remain  above  that 
threshold  for  a  minimum  of  five  sample  periods  (a  pulse  length  greater 
than  169  /xs  for  the  sampling  rate  of  31.2  IdHz) .  The  first  criteria 
represent  a  very  high  trigger  threshold,  which  was  necessary  to  avoid 
false  counts  from  sources  other  than  lightning.  Since  the  duration  of 
a  typical  lightning  flash  is  approximately  0.5  s,  the  counter  could 
have  been  triggered  several  times  by  multiple  strokes  within  a  single 
flash.  However  because  the  trigger  threshold  was  so  high,  only  a  very 

*Longltides  are  positive  east  of  Greenwich  and  negative  west, 
latitudes  are  positive  north  of  the  equator  and  negative  south. 
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small  fraction  of  the  total  strokes  triggered  the  counter.  Since  the 
first  stroke  in  a  flash  is  usually  by  far  the  strongest,  when  the 
sensor  was  triggered  it  was  by  a  first  stroke  in  most  cases.  Turman 
and  Edgar  found  that  in  approximately  10  to  20  percent  of  the  records, 
subsequent  strokes  had  also  triggered  the  counter.  They  decided  that 
since  less  than  20  percent  of  the  data  presents  any  ambiguity  between 
stroke  counting  and  flash  counting,  they  would  ignore  the  distinction 
and  assume  that  the  sensor  count  rate  was  linearly  related  to  the 
lightning  flash  rate . 

To  determine  how  many  flashes  were  detected  by  the  sensors,  the 
data  were  correlated  with  published  ground-based  records  at  several 
locations.  Various  analyses  yielded  approximately  the  same  correla¬ 
tion  factor:  the  DMSP  sensors  detected  2  percent  (±  1  percent)  of  all 
flashes  (both  I-C  and  C-G) .  Therefore  to  obtain  the  true  flash  rate, 
we  must  multiply  DMSP  map  data  by  a  factor  between  33  and  100.  Thus 
the  bias  introduced  by  the  assumption  that  the  sensor  was  triggered 
only  by  flashes  is  smaller  than  the  error  in  the  correlation  factor. 
Because  of  the  large  spread  of  these  data,  they  cannot  completely 
replace  the  old  WGL  maps;  however  they  can  be  used  to  improve  weak 
areas  of  the  old  maps  and  update  the  diurnal  modifiers. 

The  diurnal  variation  in  thunderstorm  location  can  be  observed  by 
comparing  the  "April,  Dawn"  map  (Fig.  4a)  with  "April,  Dusk" 

(Fig.  4b).  Note  that  storms  are  fairly  evenly  distributed  across  the 
land  and  ocean  masses  at  dawn,  but  concentrated  over  land  masses  at 
dusk.  Storms  during  northern  springtime  (i.e.,  April)  are  con¬ 
centrated  over  the  Northern  Hemisphere  at  dawn,  and  shift  to 
equatorial  regions  at  dusk.  Global  storm  distribution  during  southern 
spring  (northern  fall),  however,  occurs  primarily  in  the  Southern 
Hemisphere  at  both  dawn  and  dusk,  as  shown  in  Fig.  5.  Southern  spring 
distributions  also  reveal  considerably  less  storm  activity  at  dawn 
compared  to  dusk,  whereas  the  total  activity  shown  on  the  maps  of 
northern  springtime  is  approximately  the  same  for  both  dawn  and  dusk. 
The  January-March  maps  (Fig.  6)  show  a  different  diurnal  distribution 
pattern  than  either  Fig.  4  or  Fig.  5.  These  three  figures  illustrate 
the  wide  variation  in  diurnal  storm  distribution  patterns  as  well  as 
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Figure  5.  DMSP  measured  lightning  flash  occurrence  rate  maps 
for  November  through  December  at  dawn  and  dusk. 
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Figure  6.  DMSP  measured  lightning  flash  occurrence  rate  maps 
for  January  through  March  at  dawn  and  dusk. 


seasonal  variations.  They  emphasize  the  weakness  of  the  old  model, 
which  attempted  to  infer  diurnal  distribution  patterns  from  data 
measured  at  a  limited  number  of  locations. 

Although  still  insufficient  to  allow  complete  replacement  of  the 
old  TD  source  raapr ,  the  DMSP  data,  in  conjunction  with  that  from  OSO-2 
and  OSO-5  will  be  used  to  broaden  the  data  base  provided  by  the  fol¬ 
lowing  study. 

Ionosphere  Sounding  Satellite.  Identification  of  lightning 
strokes  using  radio  receivers  began  in  the  early  seventies.  Wave-form 
signatures  were  identified  that  allow  detection  of  lightning  strokes 
at  any  time  of  day.  This  technique  permits  complete  replacement  of 
the  TD  concept  with  actual  lightning  flash  density  data.  The  Iono¬ 
sphere  Sounding  Satellite  (ISS-b)  used  four  receivers  at  different 
frequencies  to  record  the  occurrence  of  global  lightning  from  June 
1978  to  May  1980.  Radio  Research  Laboratories,  Tokyo,  Japan,  pub¬ 
lished  that  information,  which,  supplemented  by  the  data  described 
above,  we  will  use  to  develop  the  new  GLO  maps.  The  new  maps  will 
replace  the  WGL  TD  maps  and  eliminate  the  need  for  the  algorithms  that 
(1)  convert  TDs  to  lightning  occurrence  and  (2)  account  for  the  diur¬ 
nal  variation  in  lightning  occurrence. 

Lightning  emissions  are  strongest  at  VLF,  and  the  intensity 
decreases  at  greater  frequencies.  However,  the  use  of  a  satellite 
platform  for  lightning  observations  requires  detection  frequencies 
higher  than  the  critical  frequency  (f^)  of  the  ionosphere  because 
waves  having  frequencies  below  f^,  are  reflected  back  toward  the 
earth.  Frequencies  jvist  above  f^,  are  most  desirable  because  the 
ionosphere  cuts  off  noise  from  oblique  angles  and  thus  limits  the 
field  of  view  (FOV)  to  the  area  directly  below  the  satellite.  Limit¬ 
ing  the  FOV  improves  source  location,  although  more  orbits  are  re¬ 
quired  to  obtain  complete  global  coverage  than  are  needed  with  a  wide- 
angle  receiver. 

The  ionosphere  varies  constantly,  so  must  be  evaluated  peri¬ 
odically  while  the  data  are  being  recorded.  Only  receivers  operating 
above  f,,  can  record  a  lightning  event.  Several  frequencies  around  the 
nominal  value  of  f^.  were  chosen  for  the  onboard  receivers  so  that  one 
detector  was  usually  operating  at  a  frequency  just  above  f^..  Those 
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frequencies  used  were  2,5,  5,  10,  and  25  MHz.  To  reduce  interference 
from  ground-based  transmitters,  frequencies  offset  above  and  below  10 
and  25  MHz  were  alternately  sampled. 

Kotaki  et  al.  [1981]  published  samples  of  the  data  from  ISS-b. 

The  lightning  counts  were  accumulated  in  geographic  bins  of  10  deg 
latitude  by  10  deg  longitude  in  six  time  blocks  of  4-h  each.  Those 
parameters  were  selected  to  correspond  to  the  standard  defined  by  CCIR 
322  [1963].  Obtaining  data  for  all  time  periods  in  the  equatorial 
regions  required  two  months  because  of  the  satellite  orbit  and  duty 
cycles  for  the  onboard  equipment.  Four  months  were  required  for  high- 
latitude  regions.  Seasonal  lightning  occurrence  maps  for  each  4-h 
time  block  were  derived. 

The  raw  count  data  for  each  bin  were  noirmalized  to  a  constant 
observation  time  for  all  sectors.  Values  for  adjacent  bins  were 
averaged  to  make  the  occurrence  rate  smooth  locally.  The  data  were 
then  converted  to  lightning  occurrence  rate  per  unit  time  and  area  as 
sho%m  in  Table  1. 

Lightning  occurrence  rate  Isopleths  for  each  of  the  six  4-h  time 
blocks  during  all  four  seasons  were  constructed.  Contour  maps  for  the 
time  block  2200  to  0200  for  each  season  are  shown  in  Figs.  7a-d.  As 
in  the  WGL  TO  maps,  seasonal  variations  in  lightning  location  are 
revealed  in  these  maps.  The  ISS-b  information  is  inherently  more 
accurate,  however,  because  the  number  of  lightning  flashes  was 
measured  directly,  whereas  WGL  used  an  imprecise  algorithm  to  convert 
TDs  to  total  number  of  lightning  flashes.  An  equally  significant 
improvement  is  the  replacement  of  the  WGL  table  of  diurnal  distribu¬ 
tion  factors  by  actxial  flash  data. 

The  ISS-b  data  represent  by  far  the  most  accurate  study  of  global 
lightning  occurrence  published  to  date.  Radio- frequency  sensors  are 
not  obscured  by  light  contamination  or  cloud  cover,  so  the  ISS-b 
detected  a  much  higher  percentage  of  total  lightning  activity  than  the 
WISP.  We  must,  however,  correlate  the  data  to  ground-based  measure¬ 
ments  before  it  can  be  used  in  the  GLO  maps . 

The  results  of  the  OSO  and  DMSP  satellite  studies  will  broaden 
the  data  base  of  the  GLO  maps  provided  by  ISS-b.  The  new  maps  will  be 
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Table  1.  Measured  rate  of  lightning  discharges  for  September 
October,  and  November  for  2200  to  0200. 
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NOTE:  Frequency  of  lightning  dischargee  occurring  within  100  km^/c  1e 

obtained  by  multiplying  the  numerical  value  by  10"  . 
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Note:  Frequency  of  lightning  discharges  occurring  within  100  km^/s 
is  obtained  by  multiplying  contour  value  by  IC®. 


ure  7.  Global  distribution  of  occurrence  rate  of  lightning 
discharges  (200  to  0200)  (Concluded). 


compared  to  CCIR  322  and  to  the  WGL  maps  to  incorporate  long-term 
variations  in  the  occurrence  rate  of  thunderstorms.  Lightning  counter 
data  originally  gathered  by  the  WGL  also  could  be  used  to  supplement 
the  GLO  map  data  base.  (That  information,  which  WGL  converted  into  TD 
form  for  their  model,  may  be  available.)  Further  studies  should  be 
incorporated  to  expand  the  time  duration  of  the  new  data  base  so  that 
variations  in  lightning  quantity  and  distribution  from  year  to  year 
can  be  included.  Such  variations  can  be  significant,  as  demonstrated 
by  the  annual  data  from  ISS-b  taken  during  the  northern  fall  season. 
Lightning  counts  for  1978  were  approximately  double  the  values  re¬ 
corded  for  1979. 

Ongoing  Measurement  Programs.  Several  organizations  have  begun 
installing  large  networks  of  lightning  counters  across  the  United 
States  to  record  the  time  and  location  of  earth  flashes.  The  most 
extensive  of  these  networks  is  operated  by  the  State  University  of  New 
York  at  Albany.  It  covers  most  of  the  Continental  United  States,  and 
new  stations  continue  to  be  added.  Although  no  data  from  this  network 
has  been  published,  portions  of  the  system  have  been  in  operation  for 
more  than  a  year,  and  some  results  may  be  obtainable.  We  will  struc¬ 
ture  the  computer  program  which  generates  the  GLO  contour  maps  to 
incorporate  new  lightning  data  as  it  becomes  available.  This  feature 
will  allow  periodic  revision  of  the  new  source  maps  by  enlarging  the 
data  base.  Ultimately,  long-term  variations  in  the  storm  occurrence 
rate  can  be  analyzed  and  incorporated. 

SOURCE  MODIFIERS. 

Noise  source  maps  describe  the  total  number  of  lightning  dis¬ 
charges  in  a  geographic  sector  over  some  specified  time  period.  As 
noted  previously  in  this  section,  different  types  of  lightning  flashes 
emit  different  frequency  spectra.  Source  modifiers  must  be  used  to 
determine  the  proportion  of  the  total  number  of  discharges  which  are 
either  I-C  or  C-G,  Although  of  lesser  significance,  sub-types  (such 
as  intercloud  and  cloud-to-air  I-C  discharges,  or  C-G  flashes  with 
positively  charged  leaders)  could  also  be  important.  If  so,  separate 
model  flashes  will  be  developed  and  additional  "submodifiers"  will  be 
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used  tc  discriminate  between  the  percentage  of  I-C  or  C-G  flashes  of 
each  subtype . 

Several  factors  have  been  identified  which  affect  the  percentage 
of  I-C  and  C-G  lightning  flashes.  The  first,  and  best  documented,  is 
latitude.  Other  factors  include  storm  severity,  whether  the  storm 
occurs  over  land  or  ocean,  and  terrain  elevation. 

Below,  we  discuss  each  factor  and  describe  the  ways  in  which  WGL 
accounted  for  them,  then  we  detail  the  approach  to  source  modifiers 
developed  for  the  new  model.  Finally,  we  review  a  modern  algorithm 
which  describes  the  variation  in  the  proportion  of  I-C  to  C-G  light¬ 
ning  with  latitude. 

WGL  Source  Modifiers. 

The  WGL  used  several  types  of  modifiers,  necessitated  in  part  by 
the  lack  of  appropriate  data  available  when  that  model  was  developed. 

Latitude .  The  variation  in  the  ratio  of  I-C  to  C-G  lightning  is 
thought  to  be  related  to  the  height  of  the  0  deg  temperature  iso¬ 
therm.  Studies  of  cloud  electrification  indicate  that  maximum  charge 
separation  occurs  near  this  level.  The  altitude  of  charge  separation 
affects  the  likelihood  that  a  discharge  will  contact  the  ground. 
Clouds  with  higher  altitude  charge  separations  produce  fewer  C-G 
flashes.  The  old  model  incorporates  the  following  equation,  which 
Pierce  [1968]  based  upon  the  scarce  data  available  at  that  time: 


(2) 


where  and  Ng  are  the  densities,  in  flashes  per  unit  area  and  time 
for  I-C  and  C-G  flashes  respectively,  and  A  is  latitude  in  degrees. 

Stoinn  Severity.  As  with  latitude,  the  basis  for  an  inverse 
relationship  between  storm  severity  (represented  by  WGL  in  the  number 
of  TDs  per  month)  and  the  percentage  of  C-G  flashes  was  the  decrease 
in  the  probability  of  C-G  strokes  as  the  altitude  of  the  charge 
separation  layer  increases.  The  strong  updrafts  associated  with 
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severe  storms  draw  low-altitude,  warm  air  higher  in  the  cloud,  in¬ 
creasing  the  altitude  of  the  charge  separation  layer.  The  WGL 
combined  the  effects  of  latitude  and  storm  severity  on  the  I-C/C-G 
ratio  as  follows ; 
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while  acknowledging  that  there  was  little  empirical  support  for  this 
modifier  equation.  Furthermore,  after  Eq.  (3)  was  used  to  calculate 
the  distribution  and  percentages  of  I-C  and  C-G  lightning  flashes  from 
TD  maps,  an  ad  hoc  correction  factor  was  applied  to  increase  the 
number  of  return  strokes  per  C-G  flash.  That  increase  was  needed  to 
make  the  predictions  of  the  WGL  model  agree  with  measurements  of 
energy  radiated  by  severe  storms . 

Occurrence  over  Land/Ocean  Masses.  The  WGL  reported  that  strikes 
to  ground  are  virtually  nonexistent  over  the  oceans.  They  accounted 
for  this  observation  by  simply  decreasing  the  WGL  model  noise  predic¬ 
tions  by  4  dB  over  the  ocean  and  2  dB  for  coastal  regions. 

Terrain  Elevation.  Although  data  indicate  that  C-G  lightning  is 
more  prevalent  in  mountainous  regions  than  other  areas,  the  effects  of 
terrain  elevation  were  not  included  in  the  WGL  noise  model.  We  will 
provide  a  modifier  to  account  for  this  effect. 

Source  Modifiers  for  the  New  Model. 

The  WGL  model  included  only  vertically  polarized  (TM)  noise. 
Because  the  emissions  from  C-G  return  strokes  dominate  TM  noise,  ad 
hoc  adjustments,  such  as  those  just  described,  could  be  used  by  WGL  to 
produce  agreement  with  experimental  data.  However,  we  must  also 
include  TE  noise,  the  major  contributor  to  which  appears  to  be  I-C 
lightning  because  of  the  long  horizontal  extent  of  its  channels.* 

The  total  number  of  lightning  strokes  (N^  +  Ng)  is  obtained  from 
the  GLO  maps.  Therefore,  an  overestimate  of  will  produce  a  cor- 

*A  two-dimensional  reconstruction  of  a  typical  I-C  flash, 
shown  in  Fig.  8,  demonstrates  that  horizontal  structure. 
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Figure  8.  Two-dimensional  reconstruction  of  intracloud  flas 
using  interferometric  lightning  location  system. 


responding  underestimate  of  Ng,  and  vice  versa.  Because  of  the  inter¬ 
dependence  of  N(,  and  Ng,  adjustments  such  as  those  used  by  WGL  cannot 
be  used.  The  new  model  must  accurately  portray  the  percentages  of 
each  lightning  type.  Since  dominates  TE  and  Ng  dominates  TM,  Nj,/Ng 
directly  affects  the  ratio  of  TE  to  TM  noise  predicted  by  the  model. 
The  predicted  TE/TM  ratio  will  be  compared  to  recent  noise  measure¬ 
ments  conducted  aboard  balloon  platforms  [Turtle  et  al . ,  in  press]. 

The  new  model  will  use  a  composite  source  definition  equation, 
which  is  the  product  of  individual  source  modifiers: 

N 

FT  -  "l-  ”s-  «e-  «g  ■ 
s 

where  Ml  is  the  latitude  factor,  Mg  is  the  storm  severity  factor,  Mg 
is  the  terrain  elevation  factor,  and  Mq  is  the  factor  which  accounts 
for  the  land/ocean  mass  ratio.  We  conducted  a  survey  of  the  published 
literature  to  develop  correlations  between  N^/Ng  and  each  of  the  above 
factors.  Below  we  discuss  the  equation  for  the  new  model's  latitude 
factor. 

Latitude .  Prentice  and  Mackerras  [1977]  derived  a  relationship 
between  the  percentages  of  I-C  lightning  discharges  and  earth  flashes 
from  29  studies  that  included  over  60,000  flashes.  The  data  points 
are  shown  in  Fig.  9.  The  discrimination  techniques  used  by  each 
experimenter  included  one  or  more  of  the  following:  visual  observa¬ 
tion,  electric  field  mill,  and  lightning  flash  counter.  The  large 
scatter  in  the  data  is  due  partly  to  the  biases  of  different  measure¬ 
ment  techniques,  but  may  also  indicate  genuine  differences  in  the  I-C 
to  C-G  ratio  for  different  localities  at  the  same  latitude  (such  as 
land/ocean  mass  differences).  Nevertheless,  they  assumed  that  such 
errors  are  random  and  that  no  systematic  bias  is  exhibited  by  the 
data. 

ITie  data  points  indicate  that  df/dA  is  small*  for  A  =  50  to 
60  deg.  Moreover,  Prentice  and  Mackerras  point  out  that  there  is  no 

*Here  we  use  C  =  N^/Ng. 
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Figure  9,  Ratio  of  densities  in  flashes  per  unit  area  and  time 
of  I-C  to  C-G  lightning  versus  latitude. 
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evidence  for  a  rapid  increase  in  N^^/Ng  as  X  approaches  zero;  and  it  is 
unlikely  that  there  are  rapid  changes  in  the  structure  and  dynamics  of 
thunderclouds  in  the  equatorial  regions.  That  assertion  is  supported 
by  Harris  and  Salman  [1972]  who  confirm  that  tropical  thunderstorm 
discharges  have  characteristics  similar  to  those  measured  in  more 
temperate  regions . 

In  view  of  the  above  considerations,  Prentice  and  Mackerras 
suggest  that  the  relationship  between  N^/Ng  and  A  must  satisfy  the 
following  boundary  conditions:  d^'/dA  =  0  at  A  =  0  deg  and  A  = 

60  deg.  They  selected  the  function  {"(A)  =  A  +  B  cos  3A  and,  by  least- 
squares  fitting,  found  A  =  4.1609  and  B  =  2.1605.  Although  better 
two-parameter  fits  were  found,  they  did  not  satisfy  the  boundary 
conditions.  Therefore  Prentice  and  Mackerras  proposed  the  empirical 
relationship 


r(A)  =  4.16  +  2.16  cos  3A  0  <  A  <  60  deg  .  (5) 

Equation  (5)  is  shown  as  the  solid  line  in  Fig.  9.  The  dashed  line  is 
from  Eq.  (2),  used  in  the  WGL  model. 

The  difference  between  the  old  and  new  relationships  is  greatest 
at  low  and  high  latitudes.  We  will  use  Eq.  (5)  for  0  <  A  <  60.  The 
WGL  omitted  sources  for  A  >  60,  We  make  the  qualitative  assumption 
that  the  structure  of  thunderstorms  at  A  >  60  does  not  vary  sig¬ 
nificantly  from  those  at  A  =  60.  For  A  >  60  we  will  use  r(60). 

Storm  Severity,  Land/Ocean  Mass,  Terrain  Elevation.  Comprehen¬ 
sive  studies,  such  as  those  described  above,  have  not  been  conducted 
to  quantify  the  effects  of  storm  severity,  land  versus  ocean  mass,  and 
terrain  elevation  on  Nj,/Ng.  Lightning  experts,  however,  generally 
acknowledge  a  relationship  between  each  of  these  factors  and  the 
proportion  of  cloud- to-earth  flashes.  Therefore,  we  will  develop 
separate  algorithms  that,  although  not  scientifically  rigorous, 
reflect  the  influence  of  each  factor  based  on  the  minimal  data  now 
available  and  modern  theoretical  understanding  of  these  phenomena.  We 
include  a  separate  modifier  for  each  factor  to  allow  easy  incorpora¬ 
tion  of  results  from  future  studies  that  better  quantify  these 
effects . 
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SECTION  3 
PROPAGATION  MODEL 


This  section  begins  by  giving  an  overview  of  PSR's  noise  propaga¬ 
tion  submodel  and  comparing  it  to  the  old  WGL  model.  It  then  gives 
examples  that  test  the  accuracy  of  the  model  against  full-wave 
calculations .  Virtually  all  mathematical  derivations  are  relegated  to 
the  appendix. 

OVERVIEW;  NEW  MODEL  VERSUS  WGL  MODEL. 

Our  method,  introduced  in  a  previous  report  [Warber  and  Field, 
1987]  ,  efficiently  calculates  mode  coupling  caused  by  changes  in 
ground  conductivity  under  ambient  day  propagation  conditions.  We 
refer  to  that  method,  and  to  the  extension  developed  in  the  present 
report,  as  the  approximate  model;  and  to  the  technique  used  by  WGL  as 
the  WGL  model.  Because  atmospheric  noise  can  seldom  be  specified  to 
within  a  few  decibels,  even  with  extensive  data,  we — somewhat 
arbitrarily — set  a  3-dB  tolerance  on  the  accuracy  of  the  propagation 
model . 

We  can  generalize  our  equation  for  the  transverse  component  of 
the  electromagnetic  field  as: 
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where  F^  is  a  transverse  field  component  (either  Ey  or  Hy)  of  the  nth 
earth- ionosphere  waveguide  mode,  P  is  the  radiated  power,  and 
(ae  sin  (x/ae))^/2  is  the  geometric  spreading  term,  where  ag  is  the 
radius  of  the  earth  and  x  the  distance  from  the  transmitter.  The  form 
of  Eq.  (6)  was  chosen  to  allow  us  to  compare  our  formulation  with  that 
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of  the  WGL.  Terms  which  normally  appear  in  the  field  equation  but 
which  are  not  important  for  this  comparison  have  been  collected  into 
K.  We  may  then  regard  k,  as  the  factor  which  converts  the  units  of  the 
rest  of  the  right-hand  side  of  Eq.  (6)  into  the  units  of  the  field  on 
the  right.  The  exact  form  of  k.  depends  on  the  kind  of  field  on  the 
right.  The  excitation  values  for  the  transmitter  and  receiver,  and 
Fj-Cx)  ,  respectively,  also  depend  on  the  type  of  field.  Note  that  F-j. 
is  a  function  of  distance  x.  The  mode  coupling  factor  is  Q^,  which 
expresses  the  scattering  of  energy  into  or  out  of  the  nth  mode.  The 
terras  ft(z)  and  fj.(x,z)  are  the  height  gain  functions  for  the  trans¬ 
mitter  and  receiver,  respectively;  they  are  functions  of  altitude  z 
and  the  receiver  term  is  also  a  function  of  distance.  can  be 

thought  of  as  the  sine  of  the  complex  angle  of  incidence  of  the  mode 
on  the  ionosphere.  Note  that  this  is  a  function  of  distance.  The 
attenuation  term  exp[-J  (x)  dx]  contains  an  integral  to  express  the 
fact  that  the  attenuation  rate  is  a  function  of  distance.  The  only 
difference  between  Eq.  (6)  and  Eq.  (57)  in  the  appendix  is  that  in  the 
latter  equation,  the  terms  that  specify  the  transmitter  and  the 
spreading  term  are  collected  into  a  single  constant  called  F^.  The 
expanded  formulation  facilitates  discussion  of  how  the  WGL  model  and 
the  approximate  model  treat  each  term. 

WGL  Model. 

To  calculate  VLF  noise,  WGL  determined  the  number  of  vertical 
lightning  strokes  within  a  15  by  15  deg  region  and  combined  those 
strokes  into  a  single  equivalent  noise  transmitter  at  the  center  of 
the  region.  The  240  equivalent  transmitters  found  in  this  way  form 
the  data  base  (the  polar  regions  are  not  included) .  For  each  trans¬ 
mitter,  the  propagation  model  is  used  to  calculate  the  signal  at  the 
center  of  the  other  239  regions.  The  RMS  sum  of  the  modal  fields  from 
all  transmitters  is  then  Che  noise  in  the  center  of  each  region. 

Using  the  RMS  sum  is  equivalent  to  assuming  that  earth-ionosphere 
waveguide  modes  do  not  interfere.  That  assumption  is  valid  because 
atmospheric  noise  is  treated  as  a  statistical  quantity  and  is  there¬ 
fore  an  average  of  the  contributions  from  numerous  sources  at 
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different  locations.  That  averaging  process  will  smear  out  nulls 
caused  by  mode  interference.  The  noise  at  intermediate  locations  is 
determined  by  interpolation. 

In  order  to  perform  this  large  number  of  calculations  quickly  on 
the  computers  available  in  1970,  the  WGL  model  retained  only  three 
modes,  whose  attenuation  and  excitation  values  were  precalculated  for 
two  ambient  ionospheres  (day  and  night)  and  several  ground 
conductivities.  The  model  takes  the  height  gain  factors  f{-(z)  and 
fj-(x,z)  to  be  unity,  since  it  deals  solely  with  TM  modes,  whose  height 
gains  are  approximately  unity  up  to  aircraft  altitudes.  The  UGL  model 
also  assumes  that  the  mode  coupling  term  Qn(x)  is  zero.  Westinghouse 
Georesearch  partially  justifies  that  assumption  by  pointing  out  that 
Qfi  is  small  for  frequencies  below  16  kHz  (they  considered  only  fre¬ 
quencies  between  10  and  30  kHz).  They  also  argued  chat  mode  coupling 
merely  rearranges  the  energy  between  modes,  at  least  during  the  sun¬ 
rise  and  sunset  periods,  and  thus  tends  to  cancel  out  in  the 
summation.  That  argument  is  incorrect,  however,  if  energy  is  coupled 
to  modes  other  than  the  lowest  three,  which  were  the  only  ones 
retained. 

The  WGL  model  includes  ad  hoc  corrections  to  account  for  the 
dependence  of  attenuation  rate  on  direction  of  propagation  and  solar 
zenith  angle.  Those  formulas  contain  constants  assumed  to  be  inde¬ 
pendent  of  the  mode  nx’mber,  although  even  WGL  points  out  that  those 
constants,  in  fact,  are  mode  dependent. 

New  Propagation  Submodel. 

The  noise  model  that  we  are  developing  will  also  use  the  concept 
of  equivalent  noise  transmitters,  although,  as  discussed  in  Sec.  2, 
our  data  base  will  include  the  effects  of  both  vertical  and  horizontal 
transmitters.  The  new  propagation  submodel  must  therefore  include 
both  TE  and  TM  modes  and  allow  for  conversion  between  them  (i.e.,  at 
night) .  The  new  model  also  must  predict  noise  under  disturbed  iono¬ 
spheric  conditions  as  well  as  ambient  conditions .  The  need  to  extend 
the  frequency  range  of  the  predictions  up  to  60  kHz  will  require  more 
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careful  treatment  of  the  mode  coupling  and  mode  conversion*,  because 
those  phenomena  are  stronger  at  LF  than  at  VLF. 

In  full-wave  computer  codes,  the  ionosphere  is  modeled  by  speci¬ 
fying  the  number-density  profiles  of  free  electrons  and  any  ion 
species  present.  These  full-wave  methods,  although  accurate,  would  be 
much  too  slow  for  the  many-source  model  of  interest  here;  a  faster 
method  is  needed.  Our  approximate  propagation  model  combines  precal¬ 
culated  full-wave  results  with  a  model  waveguide  that  treats  portions 
of  the  ionosphere  as  a  sharply  bounded  conducting  space .  That  ap¬ 
proximation,  which  is  the  major  one,  allows  formulation  of  analytic 
equations  for  mode  coupling  under  nonstratif ied  conditions.  The 
height  and  conductivity  of  the  sharply-bounded  ionosphere  are  chosen 
to  give  propagation  parameters  that  agree  with  the  full-wave  results. 

Warber  and  Field  [1987]  assumed  the  conductivity  and  height  of 
the  ionosphere  to  be  constant,  but  allowed  the  ground  conductivity  to 
change  within  a  region  of  length  Ax.  In  the  appendix  of  the  present 
report  we  extend  that  method  to  situations  in  which  the  height  of  the 
ionosphere  undergoes  lateral  changes,  but  the  conductivities  are 
constant.  Such  ionospheric  height  changes  occur,  for  example,  at 
boundaries  between  ambient  and  nuclear-disturbed  regions.  In  addi¬ 
tion,  the  appendix  derives  a  method  to  easily  calculate  the  change  in 
attenuation  rate  as  the  propagation  direction  changes  with  respect  to 
the  earth's  magnetic  field  at  night. 

The  mode  coupling  function  of  Eq.  (6)  is  the  double  summation  of 
a  series  of  scattering  terms; 
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Here  represents  the  kth  order  scattering  of  energy  from  the  nth 

to  the  mth  mode,  and  the  represents  the  energy  scattered  back 

*We  define  mode  coupling  as  the  transfer  of  energy  between 
modes  due  to  lateral  transitions  in  propagation  conditions,  whereas 
mode  conversion  results  from  the  earth's  magnetic  field  changing  part 
of  a  TM  field  into  a  TE  field  and  vice  versa. 
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into  the  mode  from  all  the  other  modes.  Warber  and  Field  [1987] 
discuss  the  scattering  matrix  in  detail.  The  computer  algorithm  that 
determines  keeps  track  of  the  size  of  each  scattering  order  and 

truncates  the  sum  when  the  scattering  becomes  small.  In  addition,  the 
number  of  retained  modes  is  monitored  and  increased,  if  necessary, 
during  a  coupling  calculation.  Examples  of  mode  coupling  at  transi¬ 
tions  between  nuclear-disturbed  and  ambient  ionospheres  are  given 
later  in  this  section. 

It  is  well-known  that  .modeling  of  propagation  under  ambient  night 
conditions  is  difficult,  especially  for  west-to-east  propagation. 

Even  full-wave  computer  codes  are  hard-pressed  to  give  accurate 
results  [Pappert  and  Hitney,  1988].  It  is  therefore  not  surprising 
that  we  were  unable  to  develop  a  satisfactory  fast-running  approximate 
method  for  undisturbed  nighttime  propagation.  However,  we  have 
developed  a  fast  method  to  calculate  the  change  in  attenuation  rate  as 
the  propagation  direction  changes.  That  method  uses  full-wave  results 
at  a  particular  propagation  direction  (e.g.,  east  to  west)  to  deter¬ 
mine  an  anisotropic  conducting  layer  to  represent  the  ionosphere. 

Then  we  can  use  the  fast  approximate  model  to  find  the  propagation 
parameters  for  other  directions,  thus  avoiding  full-fledged  calcula¬ 
tions  for  each  of  many  propagation  directions.  Again,  details  are 
given  in  the  appendix  and  later  in  this  section. 

To  summarize,  our  propagation  submodel  was  developed  to  take 
advantage  of  the  increase  in  computing  power  that  has  occurred  in  the 
18  years  since  the  publication  of  the  WGL  model.  It  is  still  neces¬ 
sary  to  perform  the  calculations  efficiently,  however,  because  of  the 
enormous  number  of  modes  and  paths.  Tests  against  the  full-wave  code 
show  our  method  decreases  computer  running  time  up  to  fivefold, 
depending  on  the  particular  situation. 

PROPAGATION  UNDER  NUCLEAR-DISTURBED  CONDITIONS. 

We  will  model  nuclear-disturbed  ionospheres  using  standard 
spread-debris  environments  based  on  radiation  from  fission  debris 
shining  down  on  the  VLF  reflecting  regions.  The  debris  is  assumed  to 
be  spread  over  a  wide  area  at  altitudes  around  150  km.  The  resulting 
disturbed  ionospheres  are  characterized  by  the  parameter  W,  which  is  a 
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measure  of  the  ionizing  intensity  of  the  debris,  and  is  defined  as 


W  = 


1  2  ’ 

A(1  +  t)  • 


(8) 


where  Yp  =  total  deposited  fission  yield  in  megatons, 

A  =  area  over  which  debris  is  uniformly  spread  in  square 
kilometers , 

t  =  time  after  the  burst  in  seconds. 


Strictly  speaking,  the  above  definition  applies  only  to  situa¬ 
tions  in  which  all  bursts  occur  at  t  -  0.  For  a  large  number  of 
bursts  at  different  times,  we  could  express  W  as  a  sum  of 
Yp/(1  +  t)^-2  terras  over  a  constant  area,  and  determine  an  equivalent 
value  for  any  time.  The  noise  model  must  apply  to  a  range  of  environ¬ 
ments  rather  than  a  specific  scenario;  the  parameter  W  is  a  simple  yet 
realistic  way  to  characterize  the  ionospheres. 

A  value  for  W  of  -  10~®  represents  a  severe  environment;  values 
between  10"^®  and  10“^^  characterize  moderate  environments;  and  smal¬ 
ler  values,  weakly  disturbed  environments. 

Figure  10  shows  the  ion  and  electron  density  profiles  calculated 
by  Field  and  Dore  [1975]  for  a  range  of  W  values.  Those  profiles  have 
been  used  in  our  full-wave  code  to  determine  the  propagation 
parameters  over  a  range  of  ground  conductivities  and  frequencies  (see 
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Vol,  II).  By  using  these  propagation  parameter  data  we  can  select  a 
set  of  effective  heights  (h)  and  ionosphere  conductivities  (aj^)  that 
cause  our  approximate  model  to  give  results  that  agree  with  full-wave 
results.  For  the  examples  that  follow,  we  found  that  at  30kHz  we 
could  take  -  3  x  10“^  S/ra  for  each  W  value  and  then  find  a  effec¬ 
tive  height  that  fit  the  full-wave  data  well.  These  are  given  in 
Table  2.  Note  that  the  effective  height  is  just  a  parameter  in  our 
model  and  does  not  correspond  to  a  physical  reflection  height. 

Moreover  such  use  of  mathematically  convenient  effective  heights  is 
common.  For  example,  when  exponential  ionospheric  models  (e.g.  0  = 
0.3,  h'  “  72  kra)  are  used  in  the  literature,  the  effective  height  h' 
merely  establishes  a  reference  level,  and  bears  no  relation  to  the 
physical  reflection  heights.  Such  reference  levels  are  not  unique  and 
can  vary  from  model  to  model  provided  that  certain  other  mathematical 
parameters  are  adjusted  to  accommodate  changes  in  h' .  For  our  model 
other  values  of  h  would  work,  provided  we  changed  accordingly. 

In  order  to  model  the  transition  from  ambient  to  disturbed  condi¬ 
tions,  we  must  specify  the  distance  over  which  this  transition  takes 
place.  This  transition  distance  is  somewhat  scenario  dependent.  For 
example,  a  high-altitude  detonation  throws  debris  high  into  the  iono¬ 
sphere,  so  at  early  times  the  debris  has  a  large  line  of  sight,  and 
the  ionosphere  is  affected  over  distances  of  several  megameters .  At 
later  times,  as  the  debris  lowers  and  spreads  out,  the  total  affected 
area  is  still  large,  but  the  width  of  transition  zone  shrinks,  because 
radiation  from  the  edge  of  the  debris  shines  over  shorter  distances. 
Roughly  speaking,  the  ionosphere  is  affected  laterally  out  to  dis¬ 
tances  on  the  order  of  3  times  the  debris  height,  often  assumed  to  be 
on  the  order  of  150  km  for  spread-fission  debris.  This  height  would 
therefore  correspond  to  a  transition  zone  width  on  the  order  of  half  a 
megameter.  More  severe  environments  will  affect  the  ionosphere  out  to 
greater  distances  than  weaker  ones. 

We  now  compare  the,  full -wave  and  approximate  calculations.  For 
the  full-wave  calculations,  we  assumed  a  transition  width  of  600  km, 
and  divided  the  transition  from  the  W  =  2  x  10~®  disturbed  ionosphere 
to  the  ambient  ionosphere  (which  has  a  W  of  approximately  2  x  10“^^) 
into  six  zones,  each  100-km  wide.  Later  in  this  section  we 
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demonstrate  the  effects  of  zone  width  on  mode  coupling.  For  the 
approximate  calculation  we  assume  the  height  of  the  ionosphere  varies 
smoothly  from  the  ambient,  h^,  to  the  disturbed  ionosphere,  hp;  i.e., 
we  assume  that  the  derivative  of  the  function  h(x)  is  continuous.  We 
define  the  coordinate  system  such  that  the  transition  region  starts  at 
X  =  0  and  ends  at  x  =  Ax.  In  the  examples  in  this  section,  we  used 
the  following  functional  form  for  h(x) : 


h(x) 


+  Ah 
+  Ah 


3(^)^ 

Ax 


for  X  >  Ax  , 

2(^)^j  for  0  <  X  <  Ax  , 
for  X  <  0  . 


(9) 


This  form  was  chosen  so  that  h(x)  and  dh/dx  are  continuous  for  all 
values  of  x.  Figure  11  illustrates  this  dependence. 

Figure  12  compares  the  full-wave  calculation  with  the  approximate 
one  for  propagation  from  a  disturbed  to  an  ambient  environment.  The 
assumed  ground  conductivity  Is  10“2  s/m.  We  show  the  field  strength 
normalized  to  unity  at  the  start  of  the  transition  region  (x  -  0  km) 
for  the  first  two  TM  modes.  For  the  approximate  calculation,  we 
assume  that  only  the  first  mode  exists  at  the  beginning  of  the 
transition.  That  assumption  is  consistent  with  the  full-wave  calcula¬ 
tion,  which  showed  that  the  second  TM  mode  was  on  the  order  of  60  dB 
smaller  than  the  first  at  the  start  of  the  transition  (1  Mm  from  the 
transmitter) .  The  solid  line  in  Fig.  12  shows  the  full-wave 
calculation.  The  Jumps  are  caused  by  mode  coupling  at  the  start  of 
each  region.  The  dashed  line  is  the  approximate  calculation.  At  the 
end  of  the  transition  (x  =«  600  km) ,  the  approximate  results  are  only 
about  1  dB  higher  than  the  full-wave  results,  well  within  the 
tolerable  error.  Similar  results  are  obtained  for  TE  modes. 

Figure  13  shows  results  for  a  signal  propagating  from  ambient 
conditions  to  disturbed  regions.  Contrary  to  the  prior  example,  the 
first  and  second  modes  dominate  at  the  start  of  the  transition,  with 
the  second,  due  to  its  greater  excitation,  being  6  dB  greater  than  the 
first.  The  second  mode  decreases  in  magnitude  and  becomes  very  small 
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Figure  11 


GruUiid 


.  Model  of  change  in  Ionosphere  height  from  disturbed 
to  ambient  day  conditions. 
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Figure  12.  Comparison  of  full-wave  to  approximate  calculation  for 
TM  propagation  from  disturbed  to  ambient  environment. 
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Figure  13.  Comparison  of  full -wave  to  approximate  calculation  for 
propagation  from  ambient  to  disturbed  environment. 
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by  the  end  of  the  transition,  whereas  the  strength  of  the  first  mode 
increases  somewhat.  This  increase  in  strength  is  due  largely  to  the 
increase  in  excitation  values  caused  by  the  narrower  waveguide.  At 
the  end  of  the  transition  the  approximate  model  results  are  1.5  dB 
lower  than  the  full-wave  results  for  the  first  mode;  again,  this  is 
within  the  3-dB  tolerance.  The  agreement  is  not  as  close  for  the 
second  mode,  but  this  mode  is  so  heavily  attenuated  that  the  errors 
are  of  no  importance  to  the  total  field. 

To  determine  the  effect  of  the  transition-zone  width  on  the 
approximate  calculation,  we  consider  propagation  from  the  disturbed  to 
the  ambient  regions  for  widths  from  0.1  to  1  Mm.  The  results  are  as 
expected;  The  shorter  transition  regions,  '-hich  have  a  larger  deriva¬ 
tive  of  h(x) ,  require  higher  scattering  orders  and  so  take  longer  to 
run.  In  Fig.  12,  for  example,  we  need  only  scattering  order  3,  while 
a  1-Mm  transition  region  needs  only  order  2,  and  a  0 , 1  Mm  region  needs 
order  7.  The  full-wave  calculations  shown  in  Figs.  12  and  13  use  four 
modes  in  regions  where  W  is  large  and  six  where  W  is  small. 

In  the  full-wave  calculation  there  is  a  trade-off  in  difficulty 
between  the  mode-finding  calculation  and  the  mode-coupling 
calculaLon.  If  the  change  in  conditions  at  a  boundary  occurs 
abruptly--or  is  severe — then  many  modes  must  be  found  in  order  to 
determine  the  signal  structure  correctly;  there  is  a  noticeable  error 
in  the  signal  if  too  few  modes  are  used.  If  more  steps  are  used  to 
make  the  change  at  the  step  boundaries  less  abrupt,  then  more  mode¬ 
coupling  calculations  are  needed. 

The  fact  that  only  order  2  is  needed  when  Ax  -  1  Mm  means  that 
very  little  mode  coupling  occurs  in  this  case.  The  change  in  the 
field  strength  of  the  dominant  mode  is  due  almost  entirely  to  the 
decrease  in  the  excitation  values  caused  by  the  increasing  height  of 
the  ionosphere.  Because  of  that  height  increase,  the  energy  in  the 
field  has  to  spread  out  over  a  greater  volume  and,  in  effect,  becomes 
diluted. 
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PROPAGATION  FOR  AMBIENT  NIGHT. 

Ambient-night  propagation  differs  from  ambient  day  or  nuclear- 
disturbed  propagation  in  three  major  ways; 

1 .  The  attenuation  rate  of  the  modes  is  much  lower . 

2.  An  individual  mode  contains  both  TM  and  TE  fields,  i.e.,  the 
modes  are  now  "quasi"  TM  or  TE. 

3.  The  propagation  parameters  are  dependent  on  the  angle  <f>  that 
the  direction  of  propagation  makes  with  the  earth's  magnetic 
field. 


In  our  model,  the  first  point  is  less  important  than  the  others, 
because  it  is  valid  to  assume  that  the  modes  do  not  interfere.  In  a 
full-wave  calculation  with  interference,  a  large  number  of  modes  (as 
many  as  20)  are  needed  to  account  for  the  signal  interference 
structure.  We  need  only  consider  six  to  eight  modes. 

The  approximate  model  treats  the  nighttime  ionosphere  as  an 
anisotropic  conducting  layer.  The  modal  equation  in  this  case  is: 


„ion  -2ikCh 
iiRii  e 
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-  xRf"'*  iiRr'* 


xRii°" 


„ion  -4ikCh 
II  Rx  e 


0  . 


(10) 


Equation  (10)  is  fully  derived  in  the  appendix. 

Here  r1°^  represents  the  ionosphere  reflection  coefficients  and 
Rgnd  (-he  ground  reflection  coefficients.  We  use  Sudden's  [1961] 
classic  definitions  in  which  iiRii  expresses  the  reflection  of  TM  fields 
into  TM  fields  and  xRx  represents  the  reflection  of  TE  fields  into  TE 
fields.  The  xRn  and  iiRx,  which  are  nonzero  only  for  the  ionosphere, 
represent  the  reflection  of  TM  fields  into  TE  fields  and  vice  versa. 
The  exp(-2ikCh)  terms  in  Eq.  (10)  occur  because  we  define  the  iono¬ 
spheric  reflection  coefficient  as  the  ratio  of  the  reflected-to- 
incident  waves  defined  at  the  ionospheric  boundary  (z  =»  h) ,  whereas 
the  ground  reflection  coefficients  are  defined  at  z  =  0.  Here  k  is 
the  wavenumber  (2fl:/wave length)  and  C  is  the  cosine  of  the  angle  of 
incidence . 
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The  height  gains  in  the  waveguide  can  be  expressed  in  terms  of 
the  reflection  coefficients.  In  Appendix  A,  we  show  how  this  is 
done . 

The  dependence  on  direction  of  propagation  can  be  dramatic. 

Figure  14  shows  the  attenuation  rate  for  several  modes  and  a  typical 
ambient-night  ionosphere.  We  chose  the  1st,  4th,  5th,  and  7th  least 
attenuated  modes*  from  a  data  set  created  at  30  kHz  for  west-to-east 
propagation.  The  assumed  ionosphere  was  exponential,  with 
/3  =  0.5  km~^  and  h'  =  80  km;  that  model  is  often  used  with  good 
results.'"'’^  The  geomagnetic  field  dip  angle  is  60  deg. 

Figure  14  shows  the  results  obtained  using  the  full-wave  code  and 
varying  the  propagation  azimuth  angle  over  360  deg.  The  lower  order 
modes  are  slighcly  less  attenuated  for  west-to-east  propagation  ((^  = 

90  deg)  than  for  east-to-west  (<^  =  270  deg)  ,  the  well-known  "east- 
west"  effect.  There  is  no  difference  between  south- to-north  (<^  =  0 
deg)  and  north- to-south  ((^  =  180  deg)  propagation,  but  there  is  a 
significant  difference  between  ^  =  0  deg  and  4>  =  ±  90  deg 
propagation.  Three  of  the  modes  show  attenuation  minima  at  <f>  -  0  deg 
or  180  deg;  one  shows  attenuation  maxima  at  those  angles.  The  three 
showing  minima  are  quasi-TM  at  <^  =  270  deg,  the  one  showing  maxima  is 
quasi-TE.  (Usually,  a  mode  that  is  quasi-TE  at  ^  =  ±  90  deg  becomes 
quasi-TM  at  ^  =  0  deg  or  180  deg,  and  vice  versa) . 

The  angle  4>  changes  along  all  paths  except  those  lying  along  the 
magnetic  equator,  or  directly  north-south.  Because  the  attenuation 
rare  depends  only  weakly  on  4>  (except  near  =  0  deg,  180  deg,  or  near 
the  magnetic  poles)  we  need  not  be  too  concerned  about  attenuation 
changing  along  a  path  due  to  changing  <f>. 

In  Appendix  A  we  show  that  in  order  to  specify  the  sharply- 
bounded  anisotropic  ionosphere,  we  need  a  parameter  y  in  addition  to 
and  h,  which  were  adequate  to  define  the  isotropic  ionosphere. 

That  additional  parameter  y,  is  given  by 

’'Note  that  the  choice  of  modes  was  arbitrary  and  made  simply  to  allow 
easy  plotting  on  a  single  figure. 

’''’'The  parameters  h'  and  p  are  the  classic  ones  originally  defined 
by  Wait  and  Spies  [1964];  also  see  Wait  [1970]. 
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Figure  14.  Attenuation  rate  as  function  of  azimuth  angle  for  four 
selected  modes  at  night. 


where  Y  and  U  are  the  standard  magnetic-ionic  parameters  defined  in 
the  appendix,  where  «g  is  the  gyro  frequency,  and  u  is  the  electron- 
neutral  ion  collision  frequency. 

In  our  approximate  model  we  use  the  modal  equation  (Eq.  10),  to 
determine  how  the  attenuation  changes  with  azimuth.  As  an  example, 
using  the  full-wave  results  for  mode  7  in  Fig.  14,  we  find  that 
h  -  122  km,’'  =  6  x  lO'^S/m,  and  2/  =  30,  give  the  best  agreement. 

We  can  then  solve  the  modal  equation  as  a  function  of  <j>  and  compare 
with  full-wave  results.  Figure  15  shows  that  comparison.  Although 
the  approximate  calculation  does  not  reproduce  all  details  of  the 
full-wave  results,  the  agreement  is  within  1  dB.  Figure  15  is  impor¬ 
tant  because  it  gives  an  alternative  to  performing  a  separate  full- 
wave  calculation  for  each  propagation  direction  from  every  noise 
source.  Instead,  we  need  perform  the  full-wave  nighttime  calculation 
only  once  for  some  cardinal  direction  and  then  use  results  like  those 
of  Fig.  15  to  adjust  to  other  directions. 


’'This  value  is  much  larger  than  the  altitude  at  which  VLF  reflections 
occur  in  the  real  ionosphere . 
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Fiqure  15.  Comparison  of  full-wavs  to  approximate  calculation  for 
selected  mode  at  night. 
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APPENDIX  A 

DERIVATION  OF  EQUATIONS 


INTRODUCTION . 

This  appendix  shows  the  derivation  of  our  method  for  rapidly- 
calculating  the  coupling  of  modes  across  a  region  of  changing  propaga¬ 
tion  conditions.  The  method  generalizes  work  of  Warber  and  Field 
[1987],  who  developed  a  model  to  compute  daytime  mode  coupling  that 
occurs  at  changes  in  ground  conductivity.  Here  we  extend  that  method 
to  lateral  changes  in  the  ionosphere,  particularly  in  the  ionosphere 
height.  Moreover,  we  account  for  the  nighttime  propagation  anisotropy 
caused  by  the  geomagnetic  field.  Some  of  the  results  of  Warber  and 
Field  [1987]  are  repeated  here  to  make  this  appendix  self-contained, 
but  the  reader  should  be  familiar  with  that  earlier  report. 

This  time  we  begin  by  deriving  the  approximate  equations  for  the 
general  case  of  a  lateral  change  in  any  propagation  condition.  Then 
we  restate  the  results  for  the  daytime  change  in  ground  conductivity 
which  was  the  subject  of  Warber  and  Field  [1987].  Next,  we  extend  the 
analysis  to  lateral  changes  in  the  ionosphere  height  which  can  occur, 
for  example,  at  the  transition  between  ambient  conditions  and  dis¬ 
turbed  ionospheric  conditions.  Finally,  we  consider  nighttime 
propagation. 

Nighttime  and  daytime  propagation  differ  in  several  important 
respects.  First,  of  course,  are  the  much  lower  attenuation  rates 
which  require  the  calculation  of  a  large  number  of  modes  in  order  to 
define  the  signal  structure.  This  is  of  less  importance  in  the  noise 
model  than  for  man-made  signals  because,  as  discussed  in  the  main 
text,  modes  from  lightning  strokes  do  not  interfere.  Of  greater 
importance  is  the  fact  that  at  night — unlike  the  daytime  case — we 
cannot  assume  that  the  modes  are  pure  TE  or  TM.  We  must  account  for  a 
single  mode  excitation  of  both  the  horizontal  magnetic  and  horizontal 
electric  fields.  Finally,  our  approximate  method  must  reproduce  the 
marked  dependence  that  the  propagation  parameters  have  on  the  angle 
between  the  propagation  direction  and  the  earth's  magnetic  field. 
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In  the  earlier  work  we  assumed  that  the  fields  in  the  ground 
could  be  ignored.  This  was  equivalent  to  assuming  that  the  Booker 
quartic  qg  in  the  ground  was  independent  of  the  propagation  parameter 
The  largest  error  this  assumption  made  in  the  value  of  qg  was 
about  2  percent  at  the  very  low  conductivity  of  sea  ice,  Ug  ~ 
lO"^  s/m.  Although  errors  of  such  small  degree  are  not  critical  to 
the  analysis,  we  do  not  make  the  same  assumption  here.  This  leads  to 
slightly  more  complicated  forms  for  the  critical  functions  and  K^, 
however,  as  defined  later  in  the  appendix.  When  those  differences 
occur  they  will  be  noted. 

We  assume  that  the  normal  to  the  wavefront  lies  in  the  xz-plane 
and  makes  an  angle  8  with  the  vertical.  We  call  the  angle  the  projec¬ 
tion  of  the  earth's  magnetic  field  on  the  ground  (xy)  plane  makes  with 
the  x-axis  the  azimuth  of  propagation  (4>) .  The  angle  the  geomagnetic 
field  makes  with  the  horizontal  is  called  the  geomagnetic  dip  and  is 
denoted  by  5.  The  dip  is  positive  in  the  Northern  Hemisphere  and 
negative  in  the  Southern.  Figure  16  shows  these  angles,  along  with 
the  coordinate  system.  The  angular  dependence  is  expressed  in  terms 
of  the  direction  cosines: 

I  =  cos  S  cos  4i  ,  (12) 

m  =  cos  6  sin  if>  ,  (13) 

n  =  -  sin  5  .  (1'^) 


By  using  Budden's  [1961]  renormalization  of  the  magnetic  field, 
the  field  H  in  this  appendix  is,  in  fact,  the  magnetic  field  multi¬ 
plied  by  the  impedance  of  free  space.  The  wave  impedances  are  there¬ 
fore  dimensionless  and  should  be  multiplied  by  377  to  put  them  into 
meter-kilogram-second  units.  The  above  normalization  also  removes  the 
distinction  between  the  magnetic  and  electric  fields.  We  often  com¬ 
bine  the  standard  TE  (Ey,  H^,  H^)  or  TM  field  components  (Hy,  Ej.,  Ez) 
into  a  generalized  field  called  F.  That  notation  allows  us  to  derive 
one  set  of  equations  that  is  applicable  to  both  TE  and  TM  modes. 
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MODEL. 


We  treat  the  earth  as  an  isotropic  conducting  half-space  below 
z  =  0,  characterized  by  the  index  of  refraction; 


2 

n 

g 


6 
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i 


(15) 


The  ionosphere  is  taken  to  be  a  possibly  anisotropic  half-space  above 
z  =  h.  The  ionosphere  is  characterized  by  index  of  refraction; 


2  ,  .  1 

n.  =  1  -  1  

1  W€ 
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(16) 


in  the  isotropic  case.  It  is  characterized  by  the  susceptibility 

A 

matrix  M  in  the  anisotropic  case ,  where 


2  2  2  2  2 

U  -  I  Y  -inYU  -  ImY  imYU  -  inY 

in YU  -  ImY^  -  m^Y^  -ilYU  -  ranY^ 

-imYU  -  InY^  ilYU  -  mnY^  -  n^Y^ 

(17) 

and  X,  Y,  and  U  are  the  normal  magneto-ionic  parameters  [see  for 
example  Ratcliffe  (1959)].  They,  and  some  other  symbols  used  in  the 
anisotropic  case,  are  defined  in  Table  3.  Note  that  in  the  isotropic 

A  A  A 

case,  where  Y/U  «  1,  M  (nj^  -  1)  I,  where  I  is  the  identity  matrix. 

We  assume  that  Cg  >  10"^  S/m,  but  make  no  assumption  about  the 
ionosphere  for  the  moment.  We  also  assume  that  the  region  where  the 
propagation  conditions  are  laterally  nonuniform  is  confined  to  the 
boundary  region;  0  <  x  <  Ax.  For  smaller  or  larger  values  of  x  the 
conditions  are  constant.  For  our  purposes,  any  changes  in  propagation 
conditions  are  considered  to  be  continuous. 


M  =  - 


X 


2  2 
U(U  -  Y  ) 
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Symbol 

w 

He 

1 ,  m,  n 

t 

m 

ZO 

N 

V 

X 

Y 

U 

C 


Table  3.  Symbols:  anisotropic  ionospheric  parameters. 

Definition 
27r  X  frequency  (s“^) 
earth's  magnetic  field  (A/m) 

direction  cosines  of  earth's  magnetic  field  (dimensionless) 
charge  of  electron  (1.6  x  10“^^  coul) 
mass  of  electron  (9.1  x  10“31 

electric  permittivity  of  free  space  (8.85  x  10~^Z  p/m) 
magnetic  permeability  of  free  space  (Att  x  10“^  H/m) 
characteristic  impedance  of  free  space  =  (3770) 

electron  density  (m“^) 
collision  frequency  of  electrons  (s"^) 

,  magneto-ionic  parameter  (dimensionless) 

magneto-ionic  parameter  (dimensionless) 

1  -  i—  ,  magneto-ionic  parameter  (dimensionless) 
w 

propagation  parameter,  cosine  of  complex  angle  of 
incidence  (dimensionless) 


M 


I 


1  -  C  ,  ImS  <  0,  sine  of  angle  of  incidence  (dimensionless) 


susceptibility  tensor  for  ionosphere  (dimensionless) 


azimuth  of  propagation  (measured  east  of  north)  (rad) 


geomagnetic  field  dip  angle  (positive  in  Northern 
Hemisphere)  (rad) 


NOTE:  All  units  are  in  the  International  System  of  Units  (SI). 
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We  asst.une  that  the  boundary  is  far  enough  from  the  transmitter 
that  the  total  field  is  made  up  of  a  sum  of  fields  of  individual 
waveguide  modes.  For  example,  the  electric  field  is 


(18) 


MAXWELL'S  EQUATIONS. 

Using  Sudden's  renormalization,  we  can  write  Maxwell's  equations 
for  the  anisotropic  case  as 


V  X  E  =  -ikH  , 


V  X  H  =  ik 


(l  +  m] 


(19) 


(20) 


Here  k  is  the  wave  number  and  is  27r /wave length.  For  an  isotropic 
ionosphere  in  the  TE  case  these  reduce  to 


dE 


H  =  , 

X  ik  Sz 


1  dE 

H  - - - - ^ 

z  ik  3x 


and 


.,2  13 

3z 


E  +  ^  H 
y  3x  z 


=  0 


For  the  TM  case.  Maxwell's  equations  become; 


ikn  E  = 

X 


3H 

dz  ’ 


(21) 


(22) 


(23) 


(24) 
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(25) 


LATERAL  CHANGES  IN  PROPAGATION  CONDITIONS:  ISOTROPIC  CASE. 

In  this  section  we  generalize  the  work  of  Warber  and  Field  [1987] 
so  that  it  applies  to  lateral  changes  in  ionospheric  height  as  well  as 
changes  in  ground  conductivity.  We  use  Fy,  the  y  component  of  the 
field,  in  the  region  where  conditions  are  not  changing,  i.e.,  x  <  0. 

(Note  that  Fy  =  Ey  in  the  TE  case  and  Fy  =  Hy  in  the  TM.)  The  form  of 

Fy  is  well  known;  i.e,, 

n 

Here  fn(z)  is  the  height  gain  function,  F^  represents  the 

forward-moving  wave,  and  Rj^  ei^^^  represents  the  backward-moving  wave, 
i.e.,  that  part  of  the  wave  reflected  from  the  boundary.  We  normalize 

the  height  gain  to  one  on  the  ground,  so  that  the  excitation  (F^)  and 

reflection  (R^)  factors  represent  the  forward  and  backward  field  on 
the  ground  at  x  =  0. 

In  the  region  where  conditions  are  spatially  changing,  we  need  a 
form  that  reduces  to  Eq.  (27).  Thus,  for  x  >  0, 

CO 

^  =  Z]  •  ^28) 

n-1 

Here  An(x)  represents  the  forward-moving  wave;  Bn(x)  represents  the 
backward-moving  wave;  and  f^Cz,  x)  is  the  height  gain,  which  is  de¬ 
pendent  on  X.  Here  An  and  Bn  are  arbitrary  functions  of  x.  We  place 
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restrictions  on  the  height-gain  term  and  determine  what  form  Aj^(x) 
and  Bj^(x)  must  have  to  satisfy  Maxwell's  equations. 

The  restriction  on  fj^  is  that  locally  it  must  satisfy  the  wave 
equation 

f  (z,  x)  =  0  ,  (29) 

n 

2  2  2 
where  q  =n(x,  z)-l+C. 

We  discuss  the  height-gain  term  fj^  in  more  detail  below  where  we 
show  that  there  exists  a  g^  that  is  orthogonal  to  f^  in  the  sense  that 


^  1  2  2,  - 
— ^  +  k  q  (x,  z) 

az 


I  g  (z,  x)f  (z,  x)  dz  =  A  (x)5 
I  °n  m  n  nm 

J  -00 


(30) 


The  boundary  conditions  satisfied  by  f^  determine  the  modal  equation 
that  gives  the  value  of  Cjj.  Thus,  the  information  about  the  propaga¬ 
tion  condition  is  contained  in  f^,  so  we  can  restrict  ourselves  to  the 
region  within  the  waveguide  (0  <  z  <  h) . 

Equations  (22)  and  (25)  can  be  written: 


F 

z 


JL  L. 

ik  dx 


F  . 

y 


(31) 


(So  in  the  TE  case  and  F^  =  -E^  in  the  TM  case.)  The  form  of 

F2  that  is  consistent  with  Eq.  (28)  and  that  reduces  to  the  correct 
form  for  x  <  0  is 


F  = 
z 


f  S 


n  n 


We  combine  Eqs .  (23)  and  (26)  (in  the  waveguide)  to  give: 


(32) 
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The  functions  A^Cx)  and  B^Cx)  must  satisfy  Eqs .  (31)  and  (33).  To 
find  out  what  this  requires  of  Ajj(x)  and  B|^(x)  ,  we  use  Eqs.  (28) 
and  (32)  ,  multiply  by  the  orthogonal  function  gj^  and  integrate  over  z 
from  -CO  to  +00.  Using  Eq.  (29)  gives,  after  some  algebra  (here,  a 
prime  symbol  denotes  the  differentiation  with  respect  to  x) , 


A'  +  B'  +  ikS 
n  n 


fA  -  B  1  +  K  fA  +  B  I  =  0  , 

n\_  n  nj  mn(^  m  mj 


(3^ 


m=l 


-  B'  +  ikS  fA  +  B  1  +  ^  fA  -  B  I  +  )  K  fA  *  BJ  —  = 
n  n  nf  n  nJ  f  n  nJ  mn(_  m  mJ 


ni=l 


Here , 


K  =  — 


ran  A 
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f:(z)  g^(z)  dz 
m  n 


Adding  Eqs.  (34)  and  (35)  yields: 
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Subtracting  those  equations  from  each  other  gives 


Now  we  define  A  and  B  in  terms  of  a  and  b,  which  are  functions  that 
vary  more  slowly  with  respect  to  x,  and  we  let 


A 

n 


a^(x) 

■  exp 

ys^(x) 


[-Ik  Xo  S„(x)  dx] 


and 


B 

n 


b  (X) 
n 

"■  exp 


S^(x)  dx]  . 


( 


Thus  Eqs.  (37)  and  (38)  become 
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(42) 


Equations  (39)  and  (40)  would  be  the  WKB  approxiraation  if  a^(x)  and 
bj^(x)  were  constants.  Thus,  Eqs .  (41)  and  (42)  becorae  the  second- 
order  WKB  approxiraation.  Up  to  this  point,  we  have  raade  no  approxi¬ 
mations.  In  Warber  and  Field  [1987],  direct  calculation  showed  that 
the  reflection  term  is  very  small  for  changes  in  ground  conductivity. 
Here  we  argue  that  the  same  is  true  for  any  other  changes  we  make. 
This  is.  at  any  rate,  the  standard  assumption  made.  We  therefore 
neglect  reflections,  which  means  we  can  set  bjjj  =  0  in  Eq.  (41). 

Removing  the  nth  term  from  the  summation  in  Eq.  (41)  yields: 


a'  +  a  K  + 
n  n  nn  2 


IE 

m=l 

nu^n 


S  +  S 
n  m 


Js  S 

m  n 


a  exp 
m  ^ 


K  =0 
ran 


(43) 


That  reduces  to  a  very  simple  form  if  we  let 

an(x)  =  a^(x)  exp  [-  d^j  -  (44) 

Then  Eq.  (43)  is 
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a;(x) 


run  m 


m=l 

nu^n 


where  for,  m  n, 


mn 


S  +  S 
n  m 


Js  S 

m  n 


[^0 


exp  iJ-  (K  -  K  )  dx 
'  ^  U  nn  mm 


exp  [-Ik  (S__  -  S„)  dx]  K^„ 


We  can  integrate  Eq.  (45)  directly,  then: 


(^5) 


(46) 


a  (X)  =  a  (0)  + 
n  n 

m=l 
rMn 

We  can  solve  Eq.  (47)  as  we  would  a  perturbation  or  scattering 
series.  That  is,  we  can  use  Eq.  (47)  to  find  am(x)  substitute 
this  into  the  integral,  getting: 


r 

Jo 


run  m 


(47) 


or  (x)  =  Q„(0)  + 
n  n 


t  VO)  j: 


7^(X)  dx 


m=l 

n»^n 


E 

m=l 

nu><n 


E 

i=i 


Jo  Jo 


(48) 
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We  can  repeat  this  process  an  arbitrary  number  of  times.  If  we  define 
7nn(x)  3  0,  and  assume  summation  over  repeated  symbols,  the  notation 
is  simpler: 


=  a  (0)  +  a"’(0)  +  a^CO)  +  G^^^  Q:"’'(0)  + 


n  nm 


nm 


nm 


(49) 


where 


™  J  n 


g;:'(x)  =  I  y(x)  dx  , 

nm  .1 Q  nm 


(50) 


G^^x)  = 


and,  in  general; 


(X) 

run 


Ei: 


(52) 


For  X  =  0,  and  ignoring  reflections,  we  have,  from  Eq.  (27): 


F 

y 


E 

•n=1 


n  n 


(53) 


At  X  =  0,  we  have,  from  Eq.  (28): 


vE 

n=1 


f  (z,  0)A  (0)  . 
n  n 


(54) 
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Equations  (53)  and  (54)  must  be  equal.  So,  since  from  Eq ,  (39)  and 
Eq.  (44), 


A^(0)  =  a^(0)/ys^(0)  =  a(Q)/JslW)  ,  (55) 

n  n  n  n  n 


we  see  that 


a  (0)  =  r  ys  (0)  . 

n  n  n 


(56) 


That  is  the  initial  condition  on  a. 

Putting  this  all  together  gives,  for  0  <  x  <  Z, 


Fy(X,  2)  = 


-  ''n[^  *  '’’‘P  (+  J'o  ‘^nn  “x] 


•  f^(x,  z)  exp  j^-ik  Jq  S^(x)  dxj 


(57) 


In  Eq.  (57),  we  define  the  mode-coupling  factor  as: 


and 


On 


E 

m=l 

nu^n 


nra 


(X) 


?  Sn(0) 


nn 


(58) 


(59) 


Here  represents  the  ith  order  scattering.  In  the  limit,  as 

goes  to  zero,  we  have: 
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CD 


(60) 


r^[i  .Q^(x)] 


E 

mssi 


nm 


VS^(O) 
S^iO) 


r  . 
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TE  HEIGHT  GAIN  AND  MODAL  EQUATIONS  FOR  ISOTROPIC  CASE. 

In  the  TE  case,  we  write  the  height-gain  function  as  e^(z) . 
Recall  that  it  is  a  solution  to  Eq.  (29).  The  boundary  conditions 
require  that  e^^  and 


A 

ik  dz  n 


(61) 


be  continuous  across  the  boundaries  at  z  =  0  and  z  =  h. 

The  solutions  of  Eq.  (29)  that  satisfy  these  conditions  as  well 
as  the  radiation  condition  (i.e.,  that  the  field  must  go  to  zero  at 
infinity),  and  have  e^(0)  =  1  are; 


e  = 
n 


Z. 

1 


fCZ^  -  11 

£ 


CZ.  +  1 

1 


-iCkh 
e  e 


-iq^k(z-h) 


cos  Ckz  +  i  sin  Ckz 
c» 


iq  kz 
.  S 


for  z  ^  h 


for  0  <  z  <  h  , 


for  z  <  0  , 


(62) 


(63) 


(64) 


where 


and 


Z  3  , 


Z.  •  —  . 


(65) 


(66) 


Recall  that  the  solution  to  Eq.  (29)  is  a  superposition  of  up-  and 
down-going  waves  whose  amplitudes  are  constants  to  be  determined  by 
the  boundary  conditions .  Thus ,  in  the  three  regions  there  are  six 
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constants.  Two  of  them  are  determined  by  the  radiation  condition  and 
one  by  the  normalization  condition.  In  reality,  the  boundary  condi¬ 
tions  at  z  =  0  and  z  =  h  determine  four  more  constants,  but  there  are 
only  three  to  be  determined.  That  fact  places  a  restriction  i^n  the 
value  of  C,  which  is  expressed  by  the  modal  equation: 


CZ 

CZ 

g 


+  1 


^^i  ~  ^  .-2iCkh 

CZ.  +  1  ® 

1 


(67) 


Note  that  we  can  write  Eq .  (67)  as 


iC  Z 

.  +  Z  cos  Ckh  = 

2  1 

1  +  C  Z.Z  sin  Ckh 

1 

1  gj  1 

1  gJ 

(68) 


The  IfLSt  form  is  particularly  easy  to  solve  numerically.  When  cTg  is 
very  large,  Zg  is  very  small,  so: 


iCZ.  cos  Ckh  =5  sin  Ckh  . 
1 


(69) 


To  zeroth  order  in  Z.  ,  C^^^w  n7r/kh,  while  to  first  order 

i  n  ' 


n  n 


1  - 


iZ^/kh 


(70) 


We  solve  Eq.  (68)  numerically  using  the  Newton-Raphson  method  starting 
with  C  =  Cj^(^). 


TM  HEIGHT  GAIN  AND  MODAL  EQUATIONS  FOR  ISOTROPIC  CASE. 

In  this  case,  the  height-gain  function  is  written  as  hn(z)  and  the 
boundary  conditions  require, that  h^  and 


1 

ikn 


_  d_ 
2  dz 


(71) 


be  continuous  across  the  boundaries.  This  gives: 
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h^(z,  x) 


f  C  -  Z 
_ g 

c  +  z. 

1 


-iCkh 

e 


Z 


-iq^k(z-h) 

e 


cos  Ckz  +  i  sin  Ckz 

c 

iq  kz 

=  S 


for  z 

for  0 

for  z 


^  h  , 

<  z  <  h 

^  0  , 


where  Zg  ^ 

Zi  =  qi/n?, 

and  the  modal  equation  is : 


C  +  Z  C  -  Z. 

_ S  =  _ i 

C  -  Z  C  +  Z. 
S  1 


-2iCkh 

e 


Equation  (77)  may  also  be  written: 


iC  |z 

.  +  Z 

cos  Ckh  =  ( 

2  1 
:  +  z.z 

sin  Ckh  . 

1 

1  gJ 

1 

1  gJ 

The  zeroth  order  solution  for  small  Zg  is  then 


n 


1 

1  +  iZ.kh 


and  the  first  order  is 


ORTHOGONAL  FUNCTION  gn(z) • 

In  Eqs .  (34)  and  (35)  we  eliminated  the  height-gain  function 

by  multiplying  by  g^  and  integrating  over  all  z.  In  this  section  we 

determine  the  form  of  gn(x)  the  TE  and  TM  cases. 

In  the  TE  case,  since  e^  is  a  solution  of  Eq.  (29),  we  can  write 

2  2 

a  a 

e  — ;r  e  -  e  — x 
m  ,  2  n  n  „  2 
5z  3z 

If  n  m,  then 


e  4-  k^fc^  -  C^l  e  e  =  0  .  (81) 

m  in  mj  n  m 
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where  the  superscripts  i,  0,  g  indicate  the  height  gain  function  in 
the  ionosphere,  the  waveguide,  and  the  ground,  respectively. 

The  radiation  conditions  ensure  that  e^,  -♦  0  as  z  -*•  ±co.  The  fact 
that  e^  and  3e^/9z  are  continuous  across  the  boundary  at  z  =  h,  and 
z  =  0  means  that  there  is  no  contribution  to  the  integral  at  those 
boundaries.  Thus,  for  the  TE  case,  g^Cz)  =  ej^(z)  . 

In  the  TM  case  9h^/3z  is  not  continuous  across  the  boundaries  at 
z  =  0,  z  =  h,  but  the  following  is: 


JL  A 

2  az 

n 


h 

n 


Thus  if  we  write'^ 


-  “2  ‘’n 
n 


(83) 


r  *  4  J 


h  h 
n  m 


dz  + 


j: 


h  h 
n  m 


dz  + 


ji: 


h  h 
n  m 


dz 


(84) 


Then  proceeding  as  in  Eq.  (82)  gives 


dz  =  0  if  n  m  . 


(85) 


’'This  point  was  dropped  inadvertently  from  Warber  and  Field 
[1987].  However,  the  results  for  and  in  the  TM  case 
do  include  the  correct  term. 
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RESULTS  FOR  LATERAL  CHANGES  OF  GROUND  CONDUCTIVITY. 

This  is  the  case  discussed  in  Warber  and  Field  [1987];  now  we 
present  the  results,  updated  to  include  the  term  that  corresponds  to 
integration  from  z  =  -oo  to  z  =  0 .  Note  that  many  quantities  we  use 
are  dependent  on  C  and  hence  are  mode  dependent.  In  order  to  cut  down 
on  the  number  of  subscripts,  we  often  suppress  the  mode  number  sub¬ 
script  on  those  quantities  that  already  have  a  subscript  such  as  qg, 
Zg,  etc.  We  do  this  only  in  those  cases  where  the  mode  number  is 
clear  from  the  content,  such  as  in  equations  that  depend  only  on  one 
mode.  We  do  use  the  mode  number  subscript  in  those  situations,  such 
as  equations  that  depend  on  two  modes,  where  it  would  be  ambiguous  not 
to  do  so. 

In  the  TE  case; 


n 


_  _1_  _ 

—  +[ikh  +  zj 

1  - 

r 

2  ■ 

► 

2ik 

»  O 

C 

.  4 

* 

(86) 


In  the  TM  case: 


-  2ik 


■  f 

^  1 

ikh  + 

2 

•  V. 

(87) 


In  Eqs .  (86)  and  (87)  the  first  term  inside  the  braces  is  due  to  the 
ground  integration.  And  for  n  ^  m,  we  can  write 


where  f(x)  =  qg  in  the  TE  case  and  {"(x)  =  Zg  in  the  TM. 
To  calculate  Knn(x)  we  note  that  in  the  TE  case. 


(88) 


nn 


i-  r*’ 

n  '*0 


(z)e  (z)  dz  . 
n 


(89) 
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However,  from  Eq.  (30)  we  see  that 


a;(x) 


2  e'(z)e  (z)  dz  . 
I  n  n 
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So, 


K 

nn 


n 


(9 


and  hence , 


£  ''nn  -  5  ' 


and  also 


(-  J'o  ^nn  H 
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In  the  TM  case , 
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*’—00  n  I  7  I  •'—CO 


dz 


-00  n 
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2  A^(x)K 
n  nn 
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21 2  2ikq 


(9: 


(9 


(9' 


(9: 


Now  the  second  term  is  much  smaller  than  1 ,  whereas  A^  -  h/2  and 

-  1.  Thus  we  can  ignore  it.  And  so  Eq.  (93)  is  good  for  the  TM 
case  also. 


RESULTS  FOR  LATERAL  CHANGES  IN  IONOSPHERE  HEIGHT, 

The  effect  of  changing  the  height  of  the  ionosphere  is 
considered.  Here  we  assume  that  h  is  a  continuous  function  of  x. 

From  either  modal  equation  [Eq3.,(67)  or  (77)]  it  is  clear  that  chang 
ing  h  changes  G,  We  can  determine  how  much  C  changes  by  taking  the 
derivative  of  the  modal  equations  with  respect  to  x.  This  gives,  in 
the  TE  case , 


C' 


-iCkh ' 


Z  +  Z.+ikh  ■ 
g  1 


(96) 


Recall  in  the  TE  case  that  Z  =  1/q.  If  we  assume  that  |Zg|  «  1  and 
|Zi|  «  kh,  which  is  true  if  aj_  2:  10"^  S/m  and  frequency  is  above 
10  kHz,  then  a  useful  approximate  solution  for  Eq .  (96)  is 


h 

^  ^  o 

C  =  C  -7-  exp 
oh 


(97) 


where  Cq  is  calculated  from  the  modal  equation  with  h  =  h^.  In  the  TM 
case , 


where 


C' 


ikh 


-  iCkh' 


1  1 


H  Z  E.Z. 
g  g  i  1 


(98) 


2  2 

H  -  C^'Cn^  +  1) 
g  g 


1;  etc. 


The  approximate  solution  of  Eq.  (98)  has  its  leading  term, 

C  =  Co  f  ,  (99) 
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which  is  the  same  as  Eq.  (97).  However,  the  range  over  which  this  is 
valid  is  not  very  useful  (typically  C  0.3,  frequency  >  40  kHz,  oi  > 
10"^  S/m)  and  so  we  will  not  use  it.  Since  Eq .  (30)  does  not  depend 
on  the  kind  of  change,  it  is  clear  that  the  are  the  same  as  in  the 
previous  case.  To  calculate  for  the  TE  case  we  write  Eq.  (36)  as 

K  =  ^ 
mn  A 

n 


e  (z)e'(z)  dz  , 
n  m 


(100) 


where  ej^(z)  is  defined  in  Eqs .  (62)  to  (64). 
To  compute  this,  note  that 


^  1  2  2 
+  k  q 
^n 


e  =  0  . 
n 


(101) 


If  we  take  the  derivative  with  respect  to  x,  then 


a^e' 


dz 
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If  n  ra,  we  can  combine  these  to  give 
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or 
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mm  n  m 


(104) 


Consider  the  last  integration  in  the  bracket: 


q  q'e  e  dz  . 
mm  m  n 


(105) 


Note  that  n^  does  not  depend  on  x  within  a  layer,  and  since  we  do  not 

move  across  boundaries  in  the  x-direction,  we  see  that 

q^  =  n^  -  1  +  c^  gives  q^Qm  =  and  thus  Eq.  (105)  becomes; 
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since  the  modes  are  orthogonal,  and  so  Eq.  (104)  becomes 
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Here  the  superscripts  i,  o,  and  g  indicate  the  field  calculated  in  the 
ionosphere,  waveguide  and  ground  respectively.  We  have  made  use  here 
of  the  continuity  of  e  and  de/dz  at  the  boundaries.  Now,  e'(z)  is 
also  continuous  at  both  boundaries,  and  de' /dz  are  continuous  at  z  = 

0.  So  Eq.  (107)  reduces  to 
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but  we  can  show  that 
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=  i  ^  k  qj  Z.  +  Z_^^  +  ikh  -  e  (h)  .  (110) 

C  im  im  gm  j  q  m 


So  substituting  these  in  Eq.  (108) 
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but  making  use  of  Eq.  (96)  gives 
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The  case  of  is  handled  in  the  same  way  as  the  earlier  case.  That 
is 
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To  calculate  for  Che  TM  case,  we  write  Eq.  (36)  as 
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dz 
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where  hji(z)  is  defined  in  Eqs .  (72)  to  (74).  As  before,  we  can  write 
equations  analogous  to  Eqs.  (102)  and  (103).  Thus  we  can  write: 
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Again  the  last  integral  reduces  to 
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which  we  have  shown  is  zero.  The  first  integral  is  done  by  breaking 
the  integration  into  the  three  ranges,  and  so: 
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Here  we  have  made  use  of  the  fact  that  h  and  1/n^  x  3h/3z  are  con¬ 
tinuous  at  the  boundaries  and  that  h(0)  =1.  As  it  turns  out,  l/n^  x 
dh' /dz  and  h'  are  also  continuous  at  z  =  0.  It  can  be  shown  that 
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and 
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Substituting  these  equations  into  Eq.  (117)  gives,  for  in  H  n: 


X  h  (h)h  (h)h’  .  (121) 

m  n 

To  find  we  can  use  the  same  logic  as  in  Eq.  (91).  Note  that,  in 

this  case,  the  conductivities  are  not  functions  of  x.  Results  using 
these  equations  are  presented  in  Sec.  3. 

NIGHT  PROPAGATION. 

In  this  section  we  discuss  propagation  under  nighttime 
ionospheres.  We  derive  a  method  to  quickly  find  the  directional 
dependence  of  the  attenuation,  and  we  derive  the  height  gain  func¬ 
tions  . 

A 

It  is  clear  from  the  susceptibility  matrix  M  (Eq.  17)  that  the 
effect  of  the  earth's  magnetic  field  on  propagation  is  contained  in 
the  ratio  of  the  magneto-ionic  parameters,  Y  and  U.  For  any  given 
frequency  Y  is  a  constant  in  the  range  25  to  150,  if  the  frequency  is 
60  to  10  kHz;  U,  however,  changes  with  altitude. 

During  ambient  day  or  disturbed  conditions,  very  low  frequency 
electromagnetic  waves  are  reflected  from  the  ionosphere  at  altitudes 
where  the  electron-neutral  collision  frequency  u  is  much  larger  than 
the  wave  frequency  w.  Hence: 


and  so. 


|U|  =  il  -  i  ^1 


(122) 
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This  means  that  the  effect  of  the  earth's  magnetic  field  is  small.  At 
night,  however,  the  waves  are  reflected  at  higher  altitudes,  and  y 
decreases  exponentially  with  a  scale  height  of  about  6.7  km.  At  about 
85  km,  y  becomes  on  the  order  of  to.  (The  exact  height  of  course 
depends  on  to.) 

If  we  write 


(123) 


where  tOg  is  called  the  gyro  frequency  and  is  about  9.3  x  10^  s~^,  then 


Y 

U 


to 

_ £_ 

to  -  iv 


(124) 


for  any  altitude  below  about  80  km.  If  we  use  the  standard  model  for 
collision  frequency; 


T  ^  T„ll  -z/(6.67  km)  -1 
i/  =  1.8xl0  e  s  , 


(125) 


then  y  is  on  the  order  of  tOg  around  66  km.  This  means  that  we  can 
regard  y  a*  Y/U  as  independent  of  frequency  below  80  km.  Let  us  look 
at  the  other  extreme  where 


to 

<j) 


»  1 


(126) 


Each  element  in  the  susceptibility  matrix  M  (Eq.  17)  has  a  denominator 
of  1  +  3/2  ~  3/2  _  Since  the  numerator  of  the  elements  of  M  are  proper- 
tional  to  either  or  y,  the  elements  of  M  reduce  to  either  0  or  a 
value  on  the  order  of  1,  depending  on  the  direction  cosines.  Again, 
these  are  independent  of  frequency.  Thus,  it  appears  we  can  find  a 
parameter  3/  which,  like  is  independent  of  frequency  and  charac¬ 
terizes  the  propagation  medium. 
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We  write  the  susceptibility  matrix  (Eq.  17)  as; 


X 

U 


3/^ 


1  - 

i.ny  -  Ino/" 


-in3/  imj/  -  Inl/ 

1  -  -iiy  -  mnl/ 


[-iin3/  -  Iny^  iiy  -  im'iP'  1  -  n^3/^ 


(127) 

A 

As  an  example  of  what  this  matrix  looks  like,  Fig.  17  plots  the  M  with 
the  -X/U  term  divided  out,  as  a  function  of  height.  We  assume  that 
the  direction  of  propagation  is  from  west  to  east  {<j>  =  90  deg)  ,  that 
the  magnetic  dip  angle  is  60  deg,  and  v  is  given  by  Eq.  (125). 

Note  that  below  about  50  and  above  80  km,  the  propagation  condi¬ 
tions  (at  least  with  respect  to  the  magnetic  effects)  are  constant, 
and  that  the  major  changes  in  conditions  occur  between  60  and  70  km, 

A 

We  plot  the  real  part  of  M  in  Fig.  17  since  the  imaginary  part  is 

A 

very  small,  and  we  want  to  emphasize  the  sign  of  the  elements  of  M. 
Although  we  have  used  30  kHz  to  plot  the  figures,  the  shape  of  the 
plots  is  in  fact  almost  independent  of  frequency. 

If  we  had  plotted  east  to  west  propagation  {ij>  =  270  deg)  then  the 
Mi3.,  M3-L,  ^23’  ^32  elements  would  change  sign.  This  causes  a 

change  in  the  value  of  the  propagation  constant  C,  which  satisfies  the 
modal  condition.  This  difference  in  C  is  reflected  in  a  very  notice¬ 
able  difference  in  the  attenuation  rate.  In  order  to  predict  this 
attenuation  within  our  propagation  models,  we  must  derive  an  expres¬ 
sion  for  the  modal  equation. 


MAXWELL'S  EQUATIONS. 

In  order  to  determine  the  modal  equation  for  night  propagation, 
we  need  to  determine  the  fields  within  an  anisotropic  layer.  As  in 
the  isotropic  case,  once  the  form  of  the  fields  in  the  ionosphere  is 
found  we  can  connect  them  to  the  fields  in  the  regions  below  the 
ionosphere  and  then  find  the  relationship  among  the  fields  that  is 
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required  for  a  mode  of  the  waveguide  to  exist.  This  relationship  is 
expressed,  of  course,  in  the  modal  equation.  To  do  this  we  use 
reflection  coefficients.  A  field  within  the  wave  guide  will  be 
reflected  at  both  the  top,  or  ionosphere,  and  bottom,  or  ground, 
boundaries.  The  reflection  coefficient  R  expresses  how  these  reflec¬ 
tions  occur.  The  requirements  of  a  mode  will  then  give  us  a  relation¬ 
ship  among  the  reflection  coefficients  that  will  be  a  form  of  the 
modal  equation. 

To  find  the  fields  in  an  anisotropic  layer,  we  will  follow  Budden 
[1964]  and  introduce  a  vector  v  whose  four  components  are  four  inde¬ 
pendent  fields  in  the  layer;  this  will  lead  to  a  simple  matrix  form 
for  Maxwell's  equations.  To  do  this,  note  that  if  we  eliminate 
from  Eqs .  (19)  and  (20)  we  can  write: 


where 


and 


>  1  a  - 

ik  az 


(128) 


V  = 


-E^ 


yj 


(129) 


[  -'""si 

SM32 

0 

l-tM33 

I+M33 

I+M33 

0 

0 

1 

0 

T  — 

M23M31 

-  «21 

- 

M23M32 

0 

SM23 

i  = 

I.M33 

H-M33 

1+M33 

1+Mii  - 

M13M31 

M13M32 

CM 

1 

0 

-SMi3 

I-M33 

I+M33 

I+M33  ^ 

(130) 


If  we  write  v(z)  =  Vq  exp(-ikqz)  then  Eq.  (128)  becomes 
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(T  -  ql)  v(z)  =  0  . 


(131) 


For  solutions  to  exist,  it  must  be  true  that 

del(T  -  ql)  =  0  .  (132) 

This  is  in  fact  Booker's  quartic  equation-  In  general,  there  are  four 
solutions  for  q.  Two  of  these  represent  up-going  waves  with  Im  q  <  0 
and  two  represent  down-going  waves  with  Im  q  >  0.  Since  we  assume 
that  the  ionosphere  is  a  conducting  half-space,  there  can  exist  no 
down-going  waves.  Thus  there  are  only  the  two  up-going  solutions 
which  we  denote  by  q^,  i  =  1,  2.  This  means  Eq.  (131)  has  two  eigen¬ 
vectors  which  we  represent  by  v(^): 


:(i) 


i  =  1,2  . 


(133) 


The  solution  of  Eq.  (132)  for  q^  and  Eq.  (131)  for  v(i)  is 
straightforward  and  does  not  concern  us  further.  Once  we  find  qj[  and 
v(l),  we  must  connect  the  results  to  the  other  layers.  To  do  this,  we 
use  Budden's  definition  of  the  reflection  coefficient.  We  next  con¬ 
sider  the  fields  to  be  represented  by  plane  waves.  That  is  any  field 
component  can  be  written  as  a  constant  times  exp(-ik(Sx+qz) ) .  The 
value  of  q  will  depend  on  the  layer  the  field  is  in,  but  Snell's  law 
says  that  the  fields  depend  on  x  in  the  same  way  in  each  layer,  i.e. 
exp(-ikSx).  We  will  suppress  this  terra  in  the  following.  Thus  the 
fields  here  are  really  just  height  gains  times  a  constant. 

The  following  notation  will  be  used  to  denote  the  fields  in  the 
ionosphere,  waveguide,  and  ground; 
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1)  the  four  fields  in  the  ionosphere  are:  E(j)(z)  and  hO)(z), 

j  =  1.2, 

2)  the  two  fields  in  the  waveguide  are:  e(z)  and  h(z) ,  and 

3)  the  two  fields  in  the  ground  are:  eS(z)  and  hS(z) . 

Let  U  represent  the  up-going  electric  field  in  the  wave  guide,  and 


u’l 

eii 

U 

ej. 


(134) 


Here  we  resolve  the  wave  into  a  component  that  is  perpendicular  to  the 
x-z  plane  ex^,  (recall  that  this  plane  contains  the  wave  normal),  and 
a  component  in  the  x-z  plane,  en^.  And,  similarly,  let 


(135) 


^  • 

We  define  the  reflection  coefficient  of  the  ionosphere  (R^®’^)  such 
that ,  at  z  =  h , 


ion 


U 


(136) 


and 


^ion 


pion 

iRll 


iRi"" 


(137) 


(In  the  following,  we  drop  the  superscript  ion  when  there  is  no  pos¬ 
sibility  of  confusion.)  Note 'that 
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(138) 


eii  = 


e  +  e 
X  z 


ex 


(139) 


and  that  the  total  field  is  the  stun  of  the  up-  and  down-going  fields. 
Consider  Fig.  18  which  shows  a  plane  wave  incident  on  the  ionosphere 
from  below  at  an  angle  8,  and  its  reflected  wave.  Since  we  have  a 
plane  wave,  it  is  easy  to  show  that  en  =  hy.  Thus, 


e  =  -  4^  h  =  -  4r  (iikCz) 

X  ik  az  y  ik  az  y 


(140) 


where  the  sign  in  the  exponent  depends  on  whether  we  have  an  up-going 
(-)  or  a  down- going  (+)  wave.  Thus  we  can  write 


total 

e 

X 


U  ^  D 
e  +  e 

X  X 


C(en 


eii) 


(141) 


also 


total 

e 

y 


U  .  D  U  D 
e  +  e  =  ei  +  ej.  . 

y  y 


From  Maxwell's  equations  we  know 


(142) 


(143) 


so , 


wU  r  U  ,  D  _  D 

h  =  -  Ce  h  =  Ce 

X  y  X  y 


(144) 


and 


X  XX 


e^  C  -h  e"^  C 

y  y 


-  C(ei  -  e?) 


(145) 


Also,  we  have 
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Figure  ’’3.  Plane  wave  incident  on  sharply  bounded  ionosphere 
and  ground  plane. 
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(1^6) 


.total  U  D 

=  eii  +  eii 

Now  using  Eqs .  (135)  and  (136)  we  can  write 

e?  =  xRx  ex  +  xRi.  ej/  (147) 

ej?  =  iiRii  eii  +  iiRx  ex  .  (148) 

The  boundary  conditions  are  that  must 

be  continuous  at  z  =  h.  Thus  we  get  at  z=h 


total 

e 

X 

=  Ceii(l  -  iiRii)  -  C||Rx  ex  -  ’ 

(149) 

total 

e 

y 

=  ex(l  +  xRx)  +  xRil  eil  =  ’ 

(150) 

•  total 
n 

X 

=  -  Cexd  -  xRx)  +  CxRii  eii  =  , 

(151) 

.total 

n 

y 

»  eii(l  +  iiRii)  +  (iRx  ex  =  , 

j*  y  ^  y 

(152) 

where  the  superscript  U  is  dropped,  since  both  en  and  ex  are  now  only 
up-going.  Here  and  K2  are  constants  that  reflect  the  fact  that  the 
eigenvectors  are  known  only  within  an  arbitrary  value.  We  can  solve 
this  set  of  equations  for  two  cases.  If  we  let  eii(h)  =  1  and  ex(h)  = 
0,  the  results  are 


iiR^" 

■H! 

L  y 

- 

K"  ■ 

e(2)1 

X 

L  y 

-  >] 

K”  ■ 

e(1>1 

X 

]■  .  (153) 

xR^" 

2C 

D 

rE<^> 

L  y 

h(2)  _  ,(2)  ^(1)1 

X  y  X  _ 

» 

(154) 
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where  the  E  and  H  field  components  are  determined  at  z  =  b,  and 


D  = 


H^2)ir 

X  j  L 


CH 


(1) 


y 


+  E 


(1) 


(155) 


If  we  let  eii(h)  =  0  and  e±(h)  =  1,  then 


=HI 

L  y 

L  y  X 

- 

CE^^^ 

L  y 

-H  E<^^ 

L  y  X 

] 

]} 


(156) 


,  ^  [£<2)  h(2)1  . 

D  L  X  y  X  y  J 


(157) 


Now  consider  the  reflection  coefficient  from  the  ground.  Since 
the  ground  is  isotropic,  there  can  be  no  conversion  for  ey  to  hy  and 
vice  versa,  so 


.  0 


(158) 


Reflection  from  the  ground  converts  a  downward- going  wave  into  an 
upward- going  wave,  so 


RU  =  D  .  (159) 


The  sitiiation  it  also  shown  in  Fig.  18.  Here  the  fields  in  the 
waveguide  at  z  »  0  are; 


total 

e 

X 


-Cen  (1  -  iiRii)  , 


(160) 
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total  D 


=  el  (1  +  xRx)  , 


total 

T 

X 


=  Cex  (1  -  iRx)  , 


.total  D  . 

=  eii  (1  +  iiRii)  . 


In  the  ground,  we  know  from  Eqs .  (63)  and  (74)  that 


e°  =  exp  (ikq^z)  , 


h^  =  exp  (ikq  z)  , 

y  S 


and  from  Maxwell's  equation  we  know  that 


h8  .  i  4- 


“  =  -rrT~e  =  qe 
X  ik  az  y  g  y 


eS  = 


_a_ 

ik  2  az  y 
g 


n 

g 


Thus  we  get  at  the  boundary,  z  =  0  (and  dropping  the  superscript 


Cej.(l  -  xRj.)  =  q  , 

O 


exd  +  xRx)  =  1  , 


and 


-  Cex(l  -  iiRii)  = 


'^g/^g  ’ 


(161) 

(162) 

(163) 

(164) 

(165) 

(166) 

(167) 

D) 

(168) 

(169) 

(170) 
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eii(l  +  ||R||)  =  1  . 


(171) 


This  gives 


and 


II  R§' 


;nd 


C  -  q  /n 

g  5 


C  +  q  /n 
g  g 


2  ’ 


iR? 


;nd 


C  +  q 


(172) 


(173) 


These  are  the  Fresnel  formulas.  Note  that  in  the  case  of  an  isotropic 
ionosphere,  we  would  have  equations  of  the  same  form  for 

For  a  mode  to  exist,  the  doubly  reflected  up-going  field  (repre¬ 
sented  on  the  right  of  Fig.  18)  must  be  the  same  as  the  initial  up- 
going  field  (on  the  left.)  This  restricts  the  allowed  values  of  C  and 
so  gives  the  modal  equation  for  the  anisotropic  case.  First  we  write 
the  up-going  field  at  any  location  in  the  waveguide  as: 


U(z) 


h^(z) 

r  u 
hy(0) 

y 

e^z) 

y  J 

II 

y 

e^(0) 

1  y  J 

^-ikqz 


(174) 


o 


(175) 


and  in  the  same  way  we  can  write  the  down-going  wave  at  any  altitude 
as 


D(z) 


B 

o 


(176) 
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Here  U  and  D  represent  constant  vectors,  U  is  arbitrary,  but 
o  o 

be  determined.  We  know 

U(h)  =  D(h)  , 


and 


D(0)  =  U(0)  . 


We  can  write  Eq.  (177)  as 


nion  rt  -ikCh  ^ 

R  U  e  =  D  e 

o  o 


Multiplying  both  sides  by  R®  and  using  Eq.  (178)  yields 


^gnd  ^ion  -  g-ikCh  ^  ^gnd  -  ^ikCh 
o  o 


+ikCh 

=  U  e 
o 


Thus  we  get: 


(^gnd  ^ion  ^-2ikCh 


A  A 

-  I)  U  =  0  . 
o 


This  requires 


det(R®’^'^ 


.ion 


-2ikCh 

e 


A 


-  I)  =  0 


or 


(iiRii 


_  gnd  - 2 ikCh 

iiRff  e 


1)  (xRr" 


_gnd  -2ikCh 
iRl  e 


1) 


_gnd  „gnd  „ion  nion  -4ihCh  „ 

-  iRj  iiRff  llRi  iRll  e  =0 


I  can 
o 


(177) 


(178) 


(179) 


(180) 


(181) 


(182) 


(183) 
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which  is  the  anisotropic  modal  equation.  This  equation  is  solved 
quickly  for  C,  once  the  parameters  of  the  ionosphere  are  found  from  a 
full  wave  calculation  in  one  direction.  Then  Eq.  (183)  is  used  to 
determine  C,  and  hence  the  attenuation  for  any  other  direction.  An 
example  of  this  is  given  in  Sec.  3.  [Note  that  in  the  case  of  an 
isotropic  ionosphere  substitution  of  Eqs .  (172)  and  (173),  along  with 

A  , 

the  similar  forms  for  would  give  the  TM  modal  equation, 

Eq.  (77),  written  in  slightly  different  form  multiplied  by  the  TE 
modal  equation,  Eq.  (67).] 


ANISOTROPIC  HEIGHT  GAIN  FUNCTIONS  e(z),  h(z)  . 

As  in  the  isotropic  case,  e(z)  and  h(z)  are  the  height  functions 
for  the  transverse  component  of  the  E  and  H  fields.  Their  form  is 


e(z) 


-iqf^^kz  -iqf^^kz 

e,  e  1  +  e„  e 


-iCkz 

e_  e  +  e,  e 

3  4 


2 

iCkz 


for  z  >  h  , 
for  0  <  z  <  h  , 


(184) 

(185) 


e 


+iq  kz 
S 


for  z  <  0  .  (186) 


and 


h(z) 


_iq(^^  k’’  kz 

e  ^'^i  +  h^  e  ^^i 


,  -iCkz  ,  iCkz 
h-  e  +  h,  e 

3  4 


_iqgkz 


for  z  >  h  , 

for  0  <  z  <  h 

for  z  <  0  . 


(187) 

(188) 

(189) 


where  q^^^ ,  q^^^  are  the  two  solutions  of  the  Booker  quartic  in  the 
ionosphere,  q|  =  n|  -  1  +  C^  ,  and  the  ej  and  hj ,  j  =1  to  4,  are 
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constants.  Note  that,  as  in  the  isotropic  case,  we  normalize  both  the 
TE  and  TM  height  gains  to  1  at  z  =  0. 

The  constants,  ej^,  62,  h;]^,  h2  can  be  determined  in  terms  of  the 
anisotropic  transmission  coefficients,  which  are  derivable  from  Eqs . 
(149)  to  (152).  However,  the  height  gains  in  the  ionosphere  are  not 
needed  to  determine  normal  nighttime  propagation  so  we  do  not  discuss 
them  further. 

To  determine  the  constants  in  the  waveguide,  we  note  that  since 
the  down-going  wave  is  in  the  waveguide. 


^iCkzl 


iCkz 

e 


Then  we  can  write 


h 


e 


4 


4 


iiRil  hj 

iRii  h^ 


-2iCkh  , 
e  + 

-2iCkh 
e  + 


II  Ri 

iRr 


636 


636 


-2iCkh 

-2iCkh 


(190) 


(191) 

(192) 


We  make  the  following  equations  less  complex  by  redefining  the 
reflection  coefficients  to  include  the  exp(-2iCkh)  term.  The  iono¬ 
sphere  reflection  coefficients  as  defined  by  Eq.  (136)  are  the  ratio 
of  fields  at  z  =  h.  Including  the  exponential  term  defines  the  coef¬ 
ficients  as  ratios  of  fields  at  z  =  0.  Thus  we  can  say  that  rI®’^  is 
referenced  on  the  ground.  We  denote  this  by  writing 
Rion^o)  m  e-2iCkh  ^  also  know  that  at  z  =  0,  we  have 

1^3  +  h^  =  1  ,  (193) 

63  +  e^  -  1  .  “  (194) 
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Thus  we  can  solve  Eqs .  (191)  to  (194)  to  give 


h3  =  |-  [l  xRi°''(0)  -  MRi°''(0)] 

63  =  ^  [l  4-  ||Rm°"(0)  -  xRil°"(0)] 
R 


where 


=  j^l  +  xRi°"(0)j[l  +  iiRn°''(0)j  -  iiRi°''(0)xRii°”(0) 


which  defines  the  height  gains  wherever  we  need  them. 


(195) 

(196) 


(197) 
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APPENDIX  B 
LIST  OF  ACRONYMS 


Acronym  Meaning  Page 


CCIR 

International  Radio  Consultative  Committee 

8 

C-G 

cloud- to-ground 

5 

DMSP 

Defense  Meteorological  Satellite  Program 

11 

ELF 

extremely  low  frequency 

1 

FOV 

field  of  view 

17 

GLO 

global  lightning  occurrence 

8 

I-C 

In-cloud 

5 

ISS 

Ionosphere  Sounding  Satellite 

17 

LF 

low  frequency 

1 

OSO 

Orbiting  Solar  Observatory 

10 

PSR 

Pacific-Sierra  Research  Corp. 

1 

RF 

radio  frequency 

4 

RMS 

root  mean  square 

30 

SI 

International  System  of  Units 

53 

TD 

thunderstorm  days 

1 

TE 

transverse  electric 

2 

TM 

transverse  magnetic 

4 

VLF 

very  low  frequency 

1 

WGL 

Westinghouse  Georesearch  Laboratory 

1 

WKB 

Wentzel-Kramers-Brillouin 

2 

WMO 

World  Meteorological  Organization 

8 
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APPENDIX  C 


LIST  OF  SYMBOLS 


In  this  appendix,  we  list  the  s}rnibols  used  in  the  text.  The  page 
number  below  indicates  the  first  page  on  which  the  symbol  is  used. 
Equation  numbers  are  also  given  when  a  quantity  is  defined  by  an 
equation.  Several  symbols  are  defined  in  Table  3  in  Appendix  A,  they 
are  also  indicated  below. 

Symbol  Description  Page 

A  -  area  over  which  debris  is  spread  (in  km^)  34 

a^  -  radius  of  the  earth  29 

Aj^(x)  -  represents  a  generalized  forward  going  wave  in 

the  transition  region  55 

a^Cx)  -  Amplitude  function  of  forward  wave  for,  Eq.  (39)  58 

Bj^(x)  -  represents  a  generalized  backward  going  wave  in 

the  transition  region  55 

b^(x)  -  an  amplitude  function  associated  with  the 

reflected  (backward  going)  wave,  Eq.  (40)  58 

Cjj(x)  -  propagation  parameter  -  cosine  of  the  complex 

eigenangle  for  mode.  50 

D  -  a  denominator,  Eq.  (155)  86 

Dg^  -  a  denominator  used  to  find  height  gains , 

Eq.  (197)  92 

D  -  a  downgoing  electric  field  in  free  space, 

Eq.  (135)  82 

E,  Ejj,  Ey,  E2  -  an  electric  field  arid  its  components.  50 

*  -  charge  of  the  electron,  Table  3  53 

e^(z)  -  TE  height  gain  function,  Eq.  (62)  63 
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Preceding  Page  Blank 


Symbol 


Description 


ei,  e2,  63,  64  -  constants  used  to  determine  the  height  gain 
of  Ey  in  the  general  anisotropic  case, 

Eq.  tl84)-(186)  90 

F  -  a  generalized  field  combining  either  the  TE 

(Ey,  Hjj,  H2)  or  TM  field  components  (Hy,  Ej^,  E^)  50 

ff,  -  critical  frequency  (MHz)  17 

fn(z)  -  generalized  height  gain  function  55 

ft-(2)  -  height  gain  factor  of  the  transmitter  29 

fj-(:<,z)  -  height  gain  factor  for  the  receiver  29 

g  (superscript)-  denotes  a  quantity  in  the  ground 

gn(z)  “  a  function  orthogonal  to  the  generalized  height 

gain  function  f^^  56 


c(j) 

nm 


(X) 


-  represents  in  some  sense  the  energy  scattered  from 
the  to  the  mode  by  the  order  scattering 
process,  Eq.  (50)  61 


H,  Hx.  Hy,  Hz 


a  renormalized  magnetic  field  and  its  components 


50 


-  geomagnetic  field  -  in  raks  units,  Table  3 


53 


^o 


-  height  of  the  ionosphere  at  start  of 
transition  zone 


36 


h][.  h2,  h3,  h4 

h(x) 

hn(z) 


constants  used  to  find  height  gain  of  Hy 

in  the  general  anisotropic  case,  Eqs .  (187)-(189) 

height  of  the  ionosphere  at  position  x 

TM  height  gain  function,  Eq.  (72) 


I  -  the  identity  matrix 

i  (superscript)-  denotes  an  ionospheric  quantity 

k  -  the  wave  numbter  =  2 n^/wave length 

Knm^^^  ~  ®  function  that  expresses  the  amount  of 

mode  coupling,  Eq.  (36) 


90 

36 

65 

52 

54 

57 
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Symbol  Description  Page 

I  -  a  direction  cosine  of  the  geomagnetic  field, 

Eq.(12)  50 

M  -  the  susceptibility  matrix  -  a  3x3  matrix,  Eq.(17)  52 

m  -  mass  of  electron.  Table  3  53 

m  -  a  direction  cosine  of  the  geomagnetic  field, 

Eq.(13)  50 

m  (subscript)  -  indicates  the  mode 

Mq  -  land/ocean  modification  factor  26 

Mg  -  terrain  elevation  modification  factor  26 

Mg  -  latitude  source  modification  factor  26 

Mg  -  storm  severity  modification  factor  26 

N  -  electron  number  density  (1/m^),  Table  3  53 

n  -  wave  normal,  see  Fig.  16  51 

n  -  a  direction  cosine  of  the  geomagnetic  field, 

Eq.  (14)  50 

n  (subscript)  -  mode  number,  Eq.  (18)  54 

N(.  -  flashes  per  unit  area  and  time  for  in-cloud 

lightning  23 

Ng  -  flashes  per  unit  area  and  time  for  cloud-to- 

ground  lightning  23 

n^,  n|,  n^  -  the  square  of  the  index  of  refraction,  in  general, 

in  the  earth,  or  in  the  ionosphere  for  the 
isotropic  case.  52 

-  Number  of  lightning  discharges/km^-month  10 

Nxd  “  Number  of  thunderstorm  days/month  10 

o  (superscript)-  denotes  a  quantity  in  free  space 

P  -  radiated  power  29 

II  (subscript)  -  indicates  a  component  parallel  to  the  x-z  plane  82 


97 


Symbol  Description  Page 

X  (subscript)  -  indicates  a  component  perpendicular  to  x-z  plane  82 

Q^(x)  -  the  total  mode-coupling  factor,  Eq .  (7)  32 

Q^'^^x)  -  the  mode-coupling  scattering  factor,  Eq.  (59)  62 

Q^-^\x)  -  represents  order  scattering  of  energy  from 

to  mode  32 

Q^^'*'^^(x)  -  represents  energy  scattered  back  into  the  mode 

from  all  other  modes,  for  order  scattering  32 

q  -  the  Booker  quartic,  q^  =  n^  -  1  +  in  an 

isotropic  region  56 

R,  R^^*^  -  the  reflection  coefficient  matrix  (a  2  x  2  matrix) , 

Eq.  (136)  82 

iiRii  >  iiRi. 

xRii ,  iiRx  ”  elements  of  the  reflection  coefficient  matrix, 

Eq.  (137)  82 

Rn  -  reflection  factor  for  mode,  indicates  the  amount 

of  reflection  of  a  wave  from  a  region  of  changing 
propagation  conditions.  55 

S^(x)  -  the  complex  sine  of  the  eigenangle  for  the 

mode  29 

t  -  time  in  seconds  34 

A 

T  -  4  X  4  matrix,  which  expresses  Maxwell's  equations 

in  an  anisotropic  region,  Eq.  (130)  80 

U  -  magneto-ionic  parameter,  Table  3  53 

U  -  an  upgoing  electric  field  in  free  space,  Eq.  (134)  82 

V  -  generalized  field  4-vector,  Eq.  (128)  80 

W  -  measure  of  the  ionizing  intensity  of  nuclear 

debris,  Eq.  (8)  34 

X  -  magneto-ionic  parameter.  Table  3  53 

X  -  distance  along  the  direction  of  propagation  52 
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z 

-  altitude 

54 

zo 

-  characteristic  impedance  of  free  space 

53 

Zi,  Zg 

-  dimensionless  impedance  for  the  ionosphere 
and  ground  respectively 

63 

-  an  excitation  function,  Eq.  (44) 

59 

P 

-  the  Wait  and  Spies  inverse  scale  height 

parameter  used  in  exponential  ionosphere  models 
as  exp(^(z-h  )). 

43 

Tn 

-  value  of  the  forward  moving  wave  on  ground 
i.e.,  the  excitation  factor 

55 

Tc-  Tr 

-  excitation  factors 

29 

-  a  cross  coupling  function,  Eq.  (46) 

60 

Ah 

-  Change  in  ionosphere  height 

37 

Ax 

-  length  of  transition  zone 

37 

5 

-  geomagnetic  field  dip  angle 

50 

^nm 

-  the  delta  function,  5^  =  1  if  n  =  m,  else  =  0 

56 

-► 

V 

-  the  gradient 

V* 

-  the  divergence 

Vx 

-  the  curl 

«0 

-  electric  permittivity  of 

free  space , 

Table  3 

«r 

-  relative  permittivity  of 

the  ground. 

Eq.  (15) 
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Symbol 

Description 

Page 

r 

- 

ratio  of  in-cloud  to  cloud -to -ground  flashes, 

Eq.  (5) 

26 

e 

- 

angle  wave  normal  makes  with  the  vertical , i .e . 
angle  of  incidence 

50 

K 

- 

a  constant, 

29 

Kl,  <2 

- 

constants 

85 

\ 

- 

latitude  in  degrees  (north  latitudes  are  positive, 
south,  are  negative) 

23 

An(x) 

- 

a  normalizing  function,  Eq.(30) 

56 

^0 

- 

magnetic  permeability  of  free  space.  Table  3 

53 

- 

C2(n2  +  1)  -  1 

70 

1/ 

- 

electron  collision  frequency  1/s,  Table  3 

53 

C7i,  Og 

- 

the  conductivity  of  the  ionosphere  and  the 
ground 

52 

r 

- 

an  integration  variable,  Eq.  (48) 

4> 

- 

azimuth  of  propagation 

50 

X 

- 

integration  variable,  used  to  denote  distance  along 
the  propagation  direction 

- 

attenuation  function 

29 

W 

- 

2jr  X  frequency 

53 

‘"g 

- 

gyro  frequency  ^  9 .3  x  10^  s"^ 

77 

(l),(2)  -  denotes  quantities  associated  with  one  of  the  two 

(superscripts)  upward-going  fields  in  an  anisotropic  ionosphere 
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